From 510436646b1bdd5b8cfab70be29334bd3cc9c828 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 29 May 2020 15:03:04 +0200 Subject: Warn if standardized residuals are unlikely normal This revealed a bug in the data returned in mkinfit$data in the case of the d_3 algorithm, which also affected the residual plot - the data from the direct fitting was not returned even if this was the better method. --- NAMESPACE | 1 + NEWS.md | 2 + R/mkinfit.R | 100 ++++++++++++---------- R/summary.mkinfit.R | 5 ++ man/mkinfit.Rd | 63 -------------- test.log | 29 ++++--- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/setup_script.R | 16 +++- tests/testthat/test_FOCUS_D_UBA_expertise.R | 10 +-- tests/testthat/test_SFORB.R | 5 ++ tests/testthat/test_analytical.R | 5 ++ tests/testthat/test_error_models.R | 17 ++-- tests/testthat/test_nafta.R | 5 +- tests/testthat/test_synthetic_data_for_UBA_2014.R | 23 +++-- 14 files changed, 131 insertions(+), 152 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index f8be4fae..0330f88f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -115,6 +115,7 @@ importFrom(stats,qnorm) importFrom(stats,qt) importFrom(stats,residuals) importFrom(stats,rnorm) +importFrom(stats,shapiro.test) importFrom(stats,update) importFrom(stats,var) importFrom(utils,getFromNamespace) diff --git a/NEWS.md b/NEWS.md index 0d5c93d3..4a923f09 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,8 @@ - 'mkinfit': Make 'use_of_ff' = "max" also the default for models specified using short names like "SFO" or "FOMC" +- 'mkinfit': Run 'stats::shapiro.test()' on standardized residuals and warn if p < 0.05 + # 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/mkinfit.R b/R/mkinfit.R index f9cb3304..d7b1b7f4 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -98,7 +98,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' @param quiet Suppress printing out the current value of the negative #' log-likelihood after each improvement? #' @param atol Absolute error tolerance, passed to [deSolve::ode()]. Default -#' is 1e-8, which is lower than the default in the [deSolve::lsoda()] +#' is 1e-8, which is lower than the default in the [deSolve::lsoda()] #' function which is used per default. #' @param rtol Absolute error tolerance, passed to [deSolve::ode()]. Default #' is 1e-10, much lower than in [deSolve::lsoda()]. @@ -119,7 +119,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' least squares fitting ("OLS") is selected. If the error model is "obs", or #' "tc", the "d_3" algorithm is selected. #' -#' The algorithm "d_3" will directly minimize the negative log-likelihood +#' The algorithm "d_3" will directly minimize the negative log-likelihood #' and independently also use the three step algorithm described below. #' The fit with the higher likelihood is returned. #' @@ -150,7 +150,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' @param trace_parms Should a trace of the parameter values be listed? #' @param \dots Further arguments that will be passed on to #' [deSolve::ode()]. -#' @importFrom stats nlminb aggregate dist +#' @importFrom stats nlminb aggregate dist shapiro.test #' @return A list with "mkinfit" in the class attribute. #' @note When using the "IORE" submodel for metabolites, fitting with #' "transform_rates = TRUE" (the default) often leads to failures of the @@ -168,7 +168,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' for measurement error in analytical chemistry. *Technometrics* 37(2), 176-184. #' #' Ranke J and Meinecke S (2019) Error Models for the Kinetic Evaluation of Chemical -#' Degradation Data. *Environments* 6(12) 124 +#' Degradation Data. *Environments* 6(12) 124 #' [doi:10.3390/environments6120124](https://doi.org/10.3390/environments6120124). #' @examples #' @@ -177,60 +177,62 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' summary(fit) #' #' # One parent compound, one metabolite, both single first order. +#' # We remove zero values from FOCUS dataset D in order to avoid warnings +#' FOCUS_D <- subset(FOCUS_2006_D, value != 0) #' # Use mkinsub for convenience in model formulation. Pathway to sink included per default. #' SFO_SFO <- mkinmod( #' parent = mkinsub("SFO", "m1"), #' m1 = mkinsub("SFO")) -#' # Fit the model to the FOCUS example dataset D using defaults -#' print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, -#' solution_type = "eigen", quiet = TRUE))) -#' parms(fit) -#' endpoints(fit) +#' +#' # Fit the model quietly to the FOCUS example dataset D using defaults +#' fit <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE) +#' # Since mkin 0.9.50.3, we get a warning about non-normality of residuals, +#' # so we try an alternative error model +#' fit.tc <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc") +#' # This avoids the warning, and the likelihood ratio test confirms it is preferable +#' lrtest(fit.tc, fit) +#' # We can also allow for different variances of parent and metabolite as error model +#' fit.obs <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "obs") +#' # This also avoids the warning about non-normality, but the two-component error model +#' # has significantly higher likelihood +#' lrtest(fit.obs, fit.tc) +#' parms(fit.tc) +#' endpoints(fit.tc) +#' +#' # We can show a quick (only one replication) benchmark for this case, as we +#' # have several alternative solution methods for the model. We skip +#' # uncompiled deSolve, as it is so slow. More benchmarks are found in the +#' # benchmark vignette #' \dontrun{ -#' # deSolve is slower when no C compiler (gcc) was available during model generation -#' print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, -#' solution_type = "deSolve"))) -#' parms(fit.deSolve) -#' endpoints(fit.deSolve) +#' if(require(rbenchmark)) { +#' benchmark(replications = 1, order = "relative", columns = c("test", "relative", "elapsed"), +#' deSolve_compiled = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc", +#' solution_type = "deSolve", use_compiled = TRUE), +#' eigen = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc", +#' solution_type = "eigen"), +#' analytical = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc", +#' solution_type = "analytical")) +#' } #' } #' -#' # Use stepwise fitting, using optimised parameters from parent only fit, FOMC +#' # Use stepwise fitting, using optimised parameters from parent only fit, FOMC-SFO #' \dontrun{ #' FOMC_SFO <- mkinmod( #' parent = mkinsub("FOMC", "m1"), #' m1 = mkinsub("SFO")) -#' # Fit the model to the FOCUS example dataset D using defaults -#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE) -#' # Use starting parameters from parent only FOMC fit -#' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE) -#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE, -#' parms.ini = fit.FOMC$bparms.ode) -#' -#' # Use stepwise fitting, using optimised parameters from parent only fit, SFORB -#' SFORB_SFO <- mkinmod( -#' parent = list(type = "SFORB", to = "m1", sink = TRUE), -#' m1 = list(type = "SFO")) -#' # Fit the model to the FOCUS example dataset D using defaults -#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, quiet = TRUE) -#' fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve", -#' quiet = TRUE) -#' # Use starting parameters from parent only SFORB fit (not really needed in this case) -#' fit.SFORB = mkinfit("SFORB", FOCUS_2006_D, quiet = TRUE) -#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, quiet = TRUE) -#' } -#' -#' \dontrun{ -#' # Weighted fits, including IRLS (error_model = "obs") -#' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), -#' m1 = mkinsub("SFO"), use_of_ff = "max") -#' f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE) -#' summary(f.noweight) -#' f.obs <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "obs", quiet = TRUE) -#' summary(f.obs) -#' f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE) -#' summary(f.tc) -#' } +#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE) +#' # Again, we get a warning and try a more sophisticated error model +#' fit.FOMC_SFO.tc <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE, error_model = "tc") +#' # This model has a higher likelihood, but not significantly so +#' lrtest(fit.tc, fit.FOMC_SFO.tc) +#' # Also, the missing standard error for log_beta and the t-tests for alpha +#' # and beta indicate overparameterisation +#' summary(fit.FOMC_SFO.tc, data = FALSE) #' +#' # We can easily use starting parameters from the parent only fit (only for illustration) +#' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE, error_model = "tc") +#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE, +#' parms.ini = fit.FOMC$bparms.ode, error_model = "tc") #' #' @export mkinfit <- function(mkinmod, observed, @@ -716,6 +718,7 @@ mkinfit <- function(mkinmod, observed, errparms <- fit_direct$par[errparms_index] } else { cost.current <- Inf # reset to avoid conflict with the OLS step + data_direct <- current_data # We need this later if it was better } } if (error_model_algorithm != "direct") { @@ -785,6 +788,7 @@ mkinfit <- function(mkinmod, observed, fit$d_3_message <- d_3_messages["direct"] degparms <- fit$par[degparms_index] errparms <- fit$par[errparms_index] + current_data <- data_direct } } } @@ -914,6 +918,10 @@ mkinfit <- function(mkinmod, observed, fit$errparms <- errparms fit$df.residual <- n_observed - length(c(degparms, errparms)) + # Check for normal distribution of residuals + fit$shapiro.p <- shapiro.test(residuals.mkinfit(fit, standardized = TRUE))$p.value + if (fit$shapiro.p < 0.05) warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05") + fit$date <- date() fit$version <- as.character(utils::packageVersion("mkin")) fit$Rversion <- paste(R.version$major, R.version$minor, sep=".") diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R index 2dc74bd7..2c291ffd 100644 --- a/R/summary.mkinfit.R +++ b/R/summary.mkinfit.R @@ -165,6 +165,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, } ans$bparms.ode <- object$bparms.ode + ans$shapiro.p <- object$shapiro.p ep <- endpoints(object) if (length(ep$ff) != 0) ans$ff <- ep$ff @@ -228,6 +229,10 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . row.names = " ")) } + if (!is.null(x$shapiro.p) && x$shapiro.p < 0.05) { + warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05") + } + cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") print(signif(x$par, digits = digits)) diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 5e6b242e..748dcb50 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -213,69 +213,6 @@ When using the "IORE" submodel for metabolites, fitting with "transform_rates = TRUE" (the default) often leads to failures of the numerical ODE solver. In this situation it may help to switch off the internal rate transformation. -} -\examples{ - -# Use shorthand notation for parent only degradation -fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) -summary(fit) - -# One parent compound, one metabolite, both single first order. -# Use mkinsub for convenience in model formulation. Pathway to sink included per default. -SFO_SFO <- mkinmod( - parent = mkinsub("SFO", "m1"), - m1 = mkinsub("SFO")) -# Fit the model to the FOCUS example dataset D using defaults -print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, - solution_type = "eigen", quiet = TRUE))) -parms(fit) -endpoints(fit) -\dontrun{ -# deSolve is slower when no C compiler (gcc) was available during model generation -print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, - solution_type = "deSolve"))) -parms(fit.deSolve) -endpoints(fit.deSolve) -} - -# Use stepwise fitting, using optimised parameters from parent only fit, FOMC -\dontrun{ -FOMC_SFO <- mkinmod( - parent = mkinsub("FOMC", "m1"), - m1 = mkinsub("SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE) -# Use starting parameters from parent only FOMC fit -fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE) -fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE, - parms.ini = fit.FOMC$bparms.ode) - -# Use stepwise fitting, using optimised parameters from parent only fit, SFORB -SFORB_SFO <- mkinmod( - parent = list(type = "SFORB", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, quiet = TRUE) -fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve", - quiet = TRUE) -# Use starting parameters from parent only SFORB fit (not really needed in this case) -fit.SFORB = mkinfit("SFORB", FOCUS_2006_D, quiet = TRUE) -fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, quiet = TRUE) -} - -\dontrun{ -# Weighted fits, including IRLS (error_model = "obs") -SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), - m1 = mkinsub("SFO"), use_of_ff = "max") -f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE) -summary(f.noweight) -f.obs <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "obs", quiet = TRUE) -summary(f.obs) -f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE) -summary(f.tc) -} - - } \references{ Rocke DM and Lorenzato S (1995) A two-component model diff --git a/test.log b/test.log index e026e559..a1fa7dfb 100644 --- a/test.log +++ b/test.log @@ -2,10 +2,10 @@ Loading mkin Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1 s] +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [0.4 s] ✔ | 7 | Fitting the SFORB model [3.3 s] -✔ | 5 | Analytical solutions for coupled models [3.1 s] +✔ | 5 | Analytical solutions for coupled models [3.2 s] ✔ | 5 | Calculation of Akaike weights ✔ | 10 | Confidence intervals and p-values [1.1 s] ✔ | 14 | Error model fitting [4.0 s] @@ -15,22 +15,29 @@ Testing mkin ✔ | 12 | Special cases of mkinfit calls [0.6 s] ✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.4 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.6 s] -✔ | 9 | Nonlinear mixed-effects models [7.7 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.5 s] +✔ | 9 | Nonlinear mixed-effects models [7.8 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] -✔ | 3 | Summary -✔ | 14 | Plotting [1.5 s] +✔ | 3 2 | Summary +──────────────────────────────────────────────────────────────────────────────── +test_plots_summary_twa.R:59: warning: Summaries are reproducible +The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05 + +test_plots_summary_twa.R:72: warning: Summaries are reproducible +The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05 +──────────────────────────────────────────────────────────────────────────────── +✔ | 14 | Plotting [1.4 s] ✔ | 4 | AIC calculation ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1 s] ✔ | 9 | Hypothesis tests [6.6 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 38.5 s +Duration: 38.3 s -OK: 159 +OK: 160 Failed: 0 -Warnings: 0 +Warnings: 2 Skipped: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index f845cc1d..7c8340cd 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2020-05-28 +Date: 2020-05-29 Optimiser: IRLS [Data] diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index def52697..8d8ba3e9 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -57,18 +57,26 @@ DFOP_SFO <- mkinmod(parent = mkinsub("DFOP", to = "m1"), # Avoid warning when fitting a dataset where zero value is removed FOCUS_D <- subset(FOCUS_2006_D, value != 0) +# We do not want warnings about non-normality of residuals here +suppressWarnings( f_sfo_sfo_desolve <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE) +) +suppressWarnings( f_sfo_sfo_eigen <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "eigen", quiet = TRUE) - +) +suppressWarnings( f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff, FOCUS_D, quiet = TRUE) +) SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data -f_2_mkin <- mkinfit("DFOP", DFOP_par_c, quiet = TRUE) +# We also suppress the warning about non-normality of residuals here, the data +# were generated with a different error model, so no wonder! +f_2_mkin <- suppressWarnings(mkinfit("DFOP", DFOP_par_c, quiet = TRUE)) f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = subset(DFOP_par_c, name == "parent")) f_2_anova <- lm(value ~ as.factor(time), data = subset(DFOP_par_c, name == "parent")) @@ -86,9 +94,9 @@ m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), fit_nw_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, quiet = TRUE) -# We know direct optimization is OK and direct needs 4 sec versus 5.5 for threestep and 6 for IRLS +# We know direct optimization is OK and direct is faster than the default d_3 fit_obs_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, error_model = "obs", quiet = TRUE, error_model_algorithm = "direct") -# We know threestep is OK, and threestep (and IRLS) need 4.8 se versus 5.6 for direct +# We know threestep is OK, and threestep (and IRLS) is faster here fit_tc_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, error_model = "tc", quiet = TRUE, error_model_algorithm = "threestep") diff --git a/tests/testthat/test_FOCUS_D_UBA_expertise.R b/tests/testthat/test_FOCUS_D_UBA_expertise.R index be2806de..101c8e15 100644 --- a/tests/testthat/test_FOCUS_D_UBA_expertise.R +++ b/tests/testthat/test_FOCUS_D_UBA_expertise.R @@ -2,11 +2,10 @@ context("Results for FOCUS D established in expertise for UBA (Ranke 2014)") # Results are from p. 40 -# Avoid warnings due to the zero value in the data for m1 at time zero -FOCUS_D <- subset(FOCUS_2006_D, value != 0) - test_that("Fits without formation fractions are correct for FOCUS D", { - fit.noff <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE) + expect_warning( + fit.noff <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE), + "p-value.*Shapiro-Wilk") expect_equal(round(as.numeric(endpoints(fit.noff)$distimes["parent", ]), 2), c(7.02, 23.33)) @@ -16,8 +15,7 @@ test_that("Fits without formation fractions are correct for FOCUS D", { }) test_that("Fits with formation fractions are correct for FOCUS D", { - skip_on_cran() - fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_D, quiet = TRUE) + fit.ff <- f_sfo_sfo.ff expect_equivalent(round(fit.ff$bparms.optim, c(2, 4, 4, 4)), c(99.60, 0.0987, 0.0053, 0.5145)) diff --git a/tests/testthat/test_SFORB.R b/tests/testthat/test_SFORB.R index 4fb736ec..91c8f2fb 100644 --- a/tests/testthat/test_SFORB.R +++ b/tests/testthat/test_SFORB.R @@ -1,5 +1,8 @@ context("Fitting the SFORB model") +# We do not want the warnings due to non-normality of residuals here +warn_option <- options(warn=-1) + test_that("Fitting the SFORB model is equivalent to fitting DFOP", { f_sforb <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE) f_dfop <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) @@ -32,3 +35,5 @@ test_that("Fitting the SFORB model is equivalent to fitting DFOP", { expect_equivalent(endpoints(f_sforb_sfo_eigen)$distimes, endpoints(f_dfop_sfo)$distimes, tolerance = 1e-6) }) + +options(warn = warn_option$warn) diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R index a34fa2cd..66fb1ace 100644 --- a/tests/testthat/test_analytical.R +++ b/tests/testthat/test_analytical.R @@ -1,5 +1,8 @@ context("Analytical solutions for coupled models") +# We do not want the warnings due to non-normality of residuals here +warn_option <- options(warn=-1) + test_that("The analytical solutions for SFO-SFO are correct", { # No sink, no formation fractions SFO_SFO_nosink <- mkinmod( @@ -58,3 +61,5 @@ test_that("The analytical solution for DFOP-SFO are correct", { tolerance = 5e-6 ) }) + +options(warn = warn_option$warn) diff --git a/tests/testthat/test_error_models.R b/tests/testthat/test_error_models.R index 2dfb2e37..169001f1 100644 --- a/tests/testthat/test_error_models.R +++ b/tests/testthat/test_error_models.R @@ -15,36 +15,37 @@ test_that("Error model 'tc' works", { test_that("The different error model fitting methods work for parent fits", { skip_on_cran() - f_9_OLS <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, - quiet = TRUE) + test_9 <- experimental_data_for_UBA_2019[[9]]$data + + f_9_OLS <- mkinfit("SFO", test_9, quiet = TRUE) expect_equivalent(round(AIC(f_9_OLS), 2), 137.43) f_9_parms_const <- parms(f_9_OLS) - f_9_direct <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_direct <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "direct", quiet = TRUE) expect_equivalent(round(AIC(f_9_direct), 2), 134.94) f_9_parms_tc_direct <- parms(f_9_direct) - f_9_twostep <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_twostep <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "twostep", quiet = TRUE) expect_equivalent(parms(f_9_twostep), f_9_parms_tc_direct, tolerance = 1e-5) - f_9_threestep <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_threestep <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "threestep", quiet = TRUE) expect_equivalent(round(AIC(f_9_threestep), 2), 139.43) expect_equivalent(parms(f_9_threestep)[1:3], f_9_parms_const) - f_9_fourstep <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_fourstep <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "fourstep", quiet = TRUE) expect_equivalent(round(AIC(f_9_fourstep), 2), 139.43) expect_equivalent(parms(f_9_fourstep)[1:3], f_9_parms_const) - f_9_IRLS <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_IRLS <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "IRLS", quiet = TRUE) expect_equivalent(round(AIC(f_9_IRLS), 2), 139.43) expect_equivalent(parms(f_9_IRLS)[1:3], f_9_parms_const) - f_9_d_3 <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + f_9_d_3 <- mkinfit("SFO", test_9, error_model = "tc", error_model_algorithm = "d_3", quiet = TRUE) expect_equivalent(round(AIC(f_9_d_3), 2), 134.94) expect_equivalent(parms(f_9_d_3), f_9_parms_tc_direct) diff --git a/tests/testthat/test_nafta.R b/tests/testthat/test_nafta.R index 2cc25f82..595db999 100644 --- a/tests/testthat/test_nafta.R +++ b/tests/testthat/test_nafta.R @@ -28,8 +28,11 @@ test_that("Test data from Appendix B are correctly evaluated", { }) test_that("Test data from Appendix D are correctly evaluated", { + # We are not interested in the warnings about non-normal residuals here + suppressWarnings( res <- nafta(NAFTA_SOP_Appendix_D, "MRID 555555", - cores = 1, quiet = TRUE) + cores = 1, quiet = TRUE) + ) # From Figure D.1 dtx_sop <- matrix(c(407, 541, 429, 1352, 5192066, 2383), nrow = 3, ncol = 2) diff --git a/tests/testthat/test_synthetic_data_for_UBA_2014.R b/tests/testthat/test_synthetic_data_for_UBA_2014.R index 4bff1b5a..989f963a 100644 --- a/tests/testthat/test_synthetic_data_for_UBA_2014.R +++ b/tests/testthat/test_synthetic_data_for_UBA_2014.R @@ -7,8 +7,8 @@ test_that("Results are correct for SFO_lin_a", { M1 = mkinsub("SFO", "M2"), M2 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) - fit_SFO_lin_a <- mkinfit(m_synth_SFO_lin, - synthetic_data_for_UBA_2014[[1]]$data, + fit_SFO_lin_a <- mkinfit(m_synth_SFO_lin, + synthetic_data_for_UBA_2014[[1]]$data, quiet = TRUE) # Results for SFO_lin_a from p. 48 @@ -21,19 +21,18 @@ test_that("Results are correct for SFO_lin_a", { # Results for DFOP_par_c from p. 54 test_that("Results are correct for DFOP_par_c", { - skip_on_cran() - m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), - M1 = mkinsub("SFO"), - M2 = mkinsub("SFO"), - use_of_ff = "max", quiet = TRUE) - + skip_on_cran() - fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par, - synthetic_data_for_UBA_2014[[12]]$data, - quiet = TRUE) + # Supress warning about non-normal residuals, the data were generated + # using a different error model, so no wonder + suppressWarnings( + fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par, + synthetic_data_for_UBA_2014[[12]]$data, + quiet = TRUE) + ) parms <- round(fit_DFOP_par_c$bparms.optim, c(1, 4, 4, 4, 4, 4, 4, 4)) - expect_equal(parms, c(parent_0 = 103.0, + expect_equal(parms, c(parent_0 = 103.0, k_M1 = 0.0389, k_M2 = 0.0095, f_parent_to_M1 = 0.5565, f_parent_to_M2 = 0.3784, k1 = 0.3263, k2 = 0.0202, g = 0.7130)) -- cgit v1.2.1