diff options
| -rw-r--r-- | R/mkinfit.R | 27 | ||||
| -rw-r--r-- | R/summary.mkinfit.R | 16 | ||||
| -rw-r--r-- | R/summary_DFOP_FOCUS_C.txt | 82 | ||||
| -rw-r--r-- | build.log | 24 | ||||
| -rw-r--r-- | man/mkinfit.Rd | 65 | ||||
| -rw-r--r-- | test.log | 31 | ||||
| -rw-r--r-- | tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt | 4 | ||||
| -rw-r--r-- | tests/testthat/summary_DFOP_FOCUS_D_eigen.txt | 4 | ||||
| -rw-r--r-- | tests/testthat/test_AIC.R | 12 | ||||
| -rw-r--r-- | tests/testthat/test_FOCUS_D_UBA_expertise.R | 2 | ||||
| -rw-r--r-- | tests/testthat/test_plot.R | 42 | ||||
| -rw-r--r-- | tests/testthat/test_plots_summary_twa.R | 128 | ||||
| -rw-r--r-- | tests/testthat/test_summary.R | 57 | ||||
| -rw-r--r-- | tests/testthat/test_twa.R | 26 |
14 files changed, 354 insertions, 166 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index d7b1b7f4..ec2d3412 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -233,7 +233,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' 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, parms.ini = "auto", @@ -258,6 +258,8 @@ mkinfit <- function(mkinmod, observed, { call <- match.call() + summary_warnings <- character() + # Derive the name used for the model if (is.character(mkinmod)) mkinmod_name <- mkinmod else mkinmod_name <- deparse(substitute(mkinmod)) @@ -289,7 +291,9 @@ mkinfit <- function(mkinmod, observed, # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model) if (any(observed$value == 0)) { - warning("Observations with value of zero were removed from the data") + zero_warning <- "Observations with value of zero were removed from the data" + summary_warnings <- c(summary_warnings, zero_warning) + warning(zero_warning) observed <- subset(observed, value != 0) } @@ -848,8 +852,9 @@ mkinfit <- function(mkinmod, observed, fit$error_model_algorithm <- error_model_algorithm if (fit$convergence != 0) { - fit$warning = paste0("Optimisation did not converge:\n", fit$message) - warning(fit$warning) + convergence_warning = paste0("Optimisation did not converge:\n", fit$message) + summary_warnings <- c(warnings, convergence_warning) + warning(convergence_warning) } else { if(!quiet) message("Optimisation successfully terminated.\n") } @@ -918,14 +923,22 @@ mkinfit <- function(mkinmod, observed, fit$errparms <- errparms fit$df.residual <- n_observed - length(c(degparms, errparms)) + # Assign the class here so method dispatch works for residuals + class(fit) <- c("mkinfit") + # 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$shapiro.p <- shapiro.test(residuals(fit, standardized = TRUE))$p.value + if (fit$shapiro.p < 0.05) { + shapiro_warning <- paste("Shapiro-Wilk test for standardized residuals: p = ", signif(fit$shapiro.p, 3)) + warning(shapiro_warning) + summary_warnings <- c(summary_warnings, shapiro_warning) + } + + fit$summary_warnings <- summary_warnings fit$date <- date() fit$version <- as.character(utils::packageVersion("mkin")) fit$Rversion <- paste(R.version$major, R.version$minor, sep=".") - class(fit) <- c("mkinfit") return(fit) } diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R index 2c291ffd..f9858c32 100644 --- a/R/summary.mkinfit.R +++ b/R/summary.mkinfit.R @@ -122,7 +122,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, date.fit = object$date, date.summary = date(), solution_type = object$solution_type, - warning = object$warning, + warnings = object$summary_warnings, use_of_ff = object$mkinmod$use_of_ff, error_model_algorithm = object$error_model_algorithm, df = c(p, rdf), @@ -190,8 +190,6 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("Date of fit: ", x$date.fit, "\n") cat("Date of summary:", x$date.summary, "\n") - if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") - cat("\nEquations:\n") nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) writeLines(strwrap(nice_diffs, exdent = 11)) @@ -223,16 +221,20 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . if(length(x$fixed$value) == 0) cat("None\n") else print(x$fixed) + # We used to only have one warning - show this for summarising old objects + if (!is.null(x[["warning"]])) cat("\n\nWarning:", x$warning, "\n\n") + + if (length(x$warnings) > 0) { + cat("\n\nWarning(s):", "\n") + cat(x$warnings, sep = "\n") + } + if (!is.null(x$AIC)) { cat("\nResults:\n\n") print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, 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/R/summary_DFOP_FOCUS_C.txt b/R/summary_DFOP_FOCUS_C.txt new file mode 100644 index 00000000..ab64a588 --- /dev/null +++ b/R/summary_DFOP_FOCUS_C.txt @@ -0,0 +1,82 @@ +mkin version used for fitting: Dummy 0.0 for testing +R version used for fitting: Dummy R version for testing +Date of fit: Dummy date for testing +Date of summary: Dummy date for testing + +Equations: +d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * + time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) + * parent + +Model predictions using solution type analytical + +Fitted using test 0 model solutions performed in test time 0 s + +Error model: Constant variance + +Error model algorithm: OLS + +Starting values for parameters to be optimised: + value type +parent_0 85.10 state +k1 0.10 deparm +k2 0.01 deparm +g 0.50 deparm + +Starting values for the transformed parameters actually optimised: + value lower upper +parent_0 85.100000 -Inf Inf +log_k1 -2.302585 -Inf Inf +log_k2 -4.605170 -Inf Inf +g_ilr 0.000000 -Inf Inf + +Fixed parameter values: +None + +Results: + + AIC BIC logLik + 29.02372 30.00984 -9.511861 + +Optimised, transformed parameters with symmetric confidence intervals: + Estimate Std. Error Lower Upper +parent_0 85.0000 0.66620 83.1500 86.8500 +log_k1 -0.7775 0.03380 -0.8713 -0.6836 +log_k2 -4.0260 0.13100 -4.3890 -3.6620 +g_ilr 1.2490 0.05811 1.0870 1.4100 +sigma 0.6962 0.16410 0.2406 1.1520 + +Parameter correlation: +[1] "Correlation matrix is platform dependent, not tested" + +Backtransformed parameters: +Confidence intervals for internally transformed parameters are asymmetric. +t-test (unrealistically) based on the assumption of normal distribution +for estimators of untransformed parameters. + Estimate t value Pr(>t) Lower Upper +parent_0 85.00000 127.600 1.131e-08 83.15000 86.85000 +k1 0.45960 29.580 3.887e-06 0.41840 0.50480 +k2 0.01785 7.636 7.901e-04 0.01241 0.02568 +g 0.85390 83.310 6.221e-08 0.82310 0.88020 +sigma 0.69620 4.243 6.618e-03 0.24060 1.15200 + +FOCUS Chi2 error levels in percent: + err.min n.optim df +All data 2.661 4 5 +parent 2.661 4 5 + +Estimated disappearance times: + DT50 DT90 DT50_k1 DT50_k2 +parent 1.887 21.25 1.508 38.83 + +Data: + time variable observed predicted residual + 0 parent 85.1 85.003 0.09726 + 1 parent 57.9 58.039 -0.13912 + 3 parent 29.9 30.054 -0.15351 + 7 parent 14.6 13.866 0.73388 + 14 parent 9.7 9.787 -0.08657 + 28 parent 6.6 7.532 -0.93205 + 63 parent 4.0 4.033 -0.03269 + 91 parent 3.9 2.447 1.45348 + 119 parent 0.6 1.484 -0.88424 @@ -2,8 +2,24 @@ * preparing ‘mkin’: * checking DESCRIPTION meta-information ... OK * installing the package to build vignettes -* creating vignettes ... OK -* checking for LF line-endings in source and make files and shell scripts -* checking for empty or unneeded directories -* building ‘mkin_0.9.50.3.tar.gz’ +* creating vignettes ... ERROR +--- re-building ‘FOCUS_D.rmd’ using rmarkdown +--- finished re-building ‘FOCUS_D.rmd’ +--- re-building ‘FOCUS_L.rmd’ using rmarkdown +Quitting from lines 80-83 (FOCUS_L.rmd) +Error: processing vignette 'FOCUS_L.rmd' failed with diagnostics: +argument 1 (type 'list') cannot be handled by 'cat' +--- failed re-building ‘FOCUS_L.rmd’ + +--- re-building ‘mkin.rmd’ using rmarkdown +--- finished re-building ‘mkin.rmd’ + +--- re-building ‘twa.rmd’ using rmarkdown +--- finished re-building ‘twa.rmd’ + +SUMMARY: processing the following file failed: + ‘FOCUS_L.rmd’ + +Error: Vignette re-building failed. +Execution halted diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 748dcb50..7f5092c5 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -214,6 +214,71 @@ When using the "IORE" submodel for metabolites, fitting with 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. +# 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 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{ +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-SFO +\dontrun{ +FOMC_SFO <- mkinmod( + parent = mkinsub("FOMC", "m1"), + m1 = mkinsub("SFO")) +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") +} +} \references{ Rocke DM and Lorenzato S (1995) A two-component model for measurement error in analytical chemistry. \emph{Technometrics} 37(2), 176-184. @@ -1,43 +1,36 @@ Loading mkin Testing mkin ✔ | OK F W S | Context +✔ | 4 | AIC calculation ✔ | 2 | Export dataset for reading into CAKE ✔ | 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] +✔ | 7 | Fitting the SFORB model [3.4 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] +✔ | 10 | Confidence intervals and p-values [1.0 s] +✔ | 14 | Error model fitting [4.1 s] ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] ✔ | 1 | Fitting the logistic model [0.2 s] ✔ | 1 | Test dataset class mkinds used in gmkin ✔ | 12 | Special cases of mkinfit calls [0.6 s] -✔ | 8 | mkinmod model generation and printing [0.2 s] +✔ | 8 | mkinmod model generation and printing [0.3 s] ✔ | 3 | Model predictions with mkinpredict [0.4 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 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 -──────────────────────────────────────────────────────────────────────────────── +✔ | 9 | Nonlinear mixed-effects models [7.7 s] ✔ | 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.5 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4 s] +✔ | 4 | Summary [0.2 s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1 s] ✔ | 9 | Hypothesis tests [6.6 s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 38.3 s +Duration: 38.2 s -OK: 160 +OK: 161 Failed: 0 -Warnings: 2 +Warnings: 0 Skipped: 0 diff --git a/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt b/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt index 9245c40b..0b31a81f 100644 --- a/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt +++ b/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt @@ -33,6 +33,10 @@ Fixed parameter values: value type m1_0 0 state + +Warning(s): +Shapiro-Wilk test for standardized residuals: p = 0.0165 + Results: AIC BIC logLik diff --git a/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt b/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt index 141f57a1..8ced5c4d 100644 --- a/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt +++ b/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt @@ -33,6 +33,10 @@ Fixed parameter values: value type m1_0 0 state + +Warning(s): +Shapiro-Wilk test for standardized residuals: p = 0.0165 + Results: AIC BIC logLik diff --git a/tests/testthat/test_AIC.R b/tests/testthat/test_AIC.R new file mode 100644 index 00000000..e9698f7c --- /dev/null +++ b/tests/testthat/test_AIC.R @@ -0,0 +1,12 @@ +context("AIC calculation") + +test_that("The AIC is reproducible", { + expect_equivalent(AIC(fits[["SFO", "FOCUS_C"]]), 59.3, scale = 1, tolerance = 0.1) + expect_equivalent(AIC(fits[, "FOCUS_C"]), + data.frame(df = c(3, 4, 5, 5), AIC = c(59.3, 44.7, 29.0, 39.2)), + scale = 1, tolerance = 0.1) + expect_error(AIC(fits["SFO", ]), "column object") + expect_equivalent(BIC(fits[, "FOCUS_C"]), + data.frame(df = c(3, 4, 5, 5), AIC = c(59.9, 45.5, 30.0, 40.2)), + scale = 1, tolerance = 0.1) +}) diff --git a/tests/testthat/test_FOCUS_D_UBA_expertise.R b/tests/testthat/test_FOCUS_D_UBA_expertise.R index 101c8e15..0a7e5219 100644 --- a/tests/testthat/test_FOCUS_D_UBA_expertise.R +++ b/tests/testthat/test_FOCUS_D_UBA_expertise.R @@ -5,7 +5,7 @@ context("Results for FOCUS D established in expertise for UBA (Ranke 2014)") test_that("Fits without formation fractions are correct for FOCUS D", { expect_warning( fit.noff <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE), - "p-value.*Shapiro-Wilk") + "Shapiro-Wilk") expect_equal(round(as.numeric(endpoints(fit.noff)$distimes["parent", ]), 2), c(7.02, 23.33)) diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R new file mode 100644 index 00000000..a33de07f --- /dev/null +++ b/tests/testthat/test_plot.R @@ -0,0 +1,42 @@ +context("Plotting") + +test_that("Plotting mkinfit and mmkin objects is reproducible", { + skip_on_cran() + plot_default_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]]) + plot_res_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]], show_residuals = TRUE) + plot_res_FOCUS_C_SFO_2 <- function() plot_res(fits[["SFO", "FOCUS_C"]]) + plot_sep_FOCUS_C_SFO <- function() plot_sep(fits[["SFO", "FOCUS_C"]]) + mkinparplot_FOCUS_C_SFO <- function() mkinparplot(fits[["SFO", "FOCUS_C"]]) + mkinerrplot_FOCUS_C_SFO <- function() mkinerrplot(fits[["SFO", "FOCUS_C"]]) + mmkin_FOCUS_C <- function() plot(fits[, "FOCUS_C"]) + mmkin_SFO <- function() plot(fits["SFO",]) + fit_D_obs_eigen <- suppressWarnings(mkinfit(SFO_SFO, FOCUS_2006_D, error_model = "obs", quiet = TRUE)) + fit_C_tc <- mkinfit("SFO", FOCUS_2006_C, error_model = "tc", quiet = TRUE) + + plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE) + plot_errmod_fit_C_tc <- function() plot_err(fit_C_tc) + + plot_res_sfo_sfo <- function() plot_res(f_sfo_sfo_desolve) + plot_err_sfo_sfo <- function() plot_err(f_sfo_sfo_desolve) + plot_errmod_fit_obs_1 <- function() plot_err(fit_obs_1, sep_obs = FALSE) + plot_errmod_fit_tc_1 <- function() plot_err(fit_tc_1, sep_obs = FALSE) + + skip_if(getRversion() > 4.0) + vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with defaults", plot_default_FOCUS_C_SFO) + vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with residuals like in gmkin", plot_res_FOCUS_C_SFO) + vdiffr::expect_doppelganger("plot_res for FOCUS C", plot_res_FOCUS_C_SFO_2) + vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with sep = TRUE", plot_sep_FOCUS_C_SFO) + vdiffr::expect_doppelganger("mkinparplot for FOCUS C SFO", mkinparplot_FOCUS_C_SFO) + vdiffr::expect_doppelganger("mkinerrplot for FOCUS C SFO", mkinerrplot_FOCUS_C_SFO) + vdiffr::expect_doppelganger("mmkin plot for FOCUS C", mmkin_FOCUS_C) + vdiffr::expect_doppelganger("mmkin plot for SFO (FOCUS C and D)", mmkin_SFO) + vdiffr::expect_doppelganger("plot_errmod with FOCUS C tc", plot_errmod_fit_C_tc) + skip_on_travis() # Still not working on Travis, maybe because of deSolve producing + # different results when not working with a compiled model or eigenvalues + vdiffr::expect_doppelganger("plot_errmod with FOCUS D obs eigen", plot_errmod_fit_D_obs_eigen) + vdiffr::expect_doppelganger("plot_res for FOCUS D", plot_res_sfo_sfo) + vdiffr::expect_doppelganger("plot_err for FOCUS D", plot_err_sfo_sfo) + vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_tc", plot_errmod_fit_tc_1) + vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_obs", plot_errmod_fit_obs_1) +}) + diff --git a/tests/testthat/test_plots_summary_twa.R b/tests/testthat/test_plots_summary_twa.R deleted file mode 100644 index aedc9da3..00000000 --- a/tests/testthat/test_plots_summary_twa.R +++ /dev/null @@ -1,128 +0,0 @@ -context("Calculation of maximum time weighted average concentrations (TWAs)") - -test_that("Time weighted average concentrations are correct", { - skip_on_cran() - - outtimes_10 <- seq(0, 10, length.out = 10000) - - ds <- "FOCUS_C" - for (model in models) { - fit <- fits[[model, ds]] - bpar <- summary(fit)$bpar[, "Estimate"] - pred_10 <- mkinpredict(fit$mkinmod, - odeparms = bpar[2:length(bpar)], - odeini = c(parent = bpar[[1]]), - outtimes = outtimes_10) - twa_num <- mean(pred_10[, "parent"]) - names(twa_num) <- 10 - twa_ana <- max_twa_parent(fit, 10) - - # Test for absolute difference (scale = 1) - # The tolerance can be reduced if the length of outtimes is increased, - # but this needs more computing time so we stay with lenght.out = 10k - expect_equal(twa_num, twa_ana, tolerance = 0.003, scale = 1) - } -}) - -context("Summary") - -test_that("Summaries are reproducible", { - fit <- fits[["DFOP", "FOCUS_C"]] - test_summary <- summary(fit) - test_summary$fit_version <- "Dummy 0.0 for testing" - test_summary$fit_Rversion <- "Dummy R version for testing" - test_summary$date.fit <- "Dummy date for testing" - test_summary$date.summary <- "Dummy date for testing" - test_summary$calls <- "test 0" - test_summary$Corr <- signif(test_summary$Corr, 1) - test_summary$time <- c(elapsed = "test time 0") - # The correlation matrix is quite platform dependent - # It differs between i386 and amd64 on Windows - # and between Travis and my own Linux system - test_summary$Corr <- "Correlation matrix is platform dependent, not tested" - expect_known_output(print(test_summary), "summary_DFOP_FOCUS_C.txt") - - test_summary_2 <- summary(f_sfo_sfo_eigen) - test_summary_2$fit_version <- "Dummy 0.0 for testing" - test_summary_2$fit_Rversion <- "Dummy R version for testing" - test_summary_2$date.fit <- "Dummy date for testing" - test_summary_2$date.summary <- "Dummy date for testing" - test_summary_2$calls <- "test 0" - test_summary_2$time <- c(elapsed = "test time 0") - # The correlation matrix is quite platform dependent - # It differs between i386 and amd64 on Windows - # and between Travis and my own Linux system - # Even more so when using the Eigen method - test_summary_2$Corr <- "Correlation matrix is platform dependent, not tested" - # The residuals for this method are also platform sensitive - test_summary_2$data$residual <- "not tested" - expect_known_output(print(test_summary_2), "summary_DFOP_FOCUS_D_eigen.txt") - - test_summary_3 <- summary(f_sfo_sfo_desolve) - test_summary_3$fit_version <- "Dummy 0.0 for testing" - test_summary_3$fit_Rversion <- "Dummy R version for testing" - test_summary_3$date.fit <- "Dummy date for testing" - test_summary_3$date.summary <- "Dummy date for testing" - test_summary_3$calls <- "test 0" - test_summary_3$time <- c(elapsed = "test time 0") - # The correlation matrix is quite platform dependent - # It differs between i386 and amd64 on Windows - # and between Travis and my own Linux system - test_summary_3$Corr <- "Correlation matrix is platform dependent, not tested" - expect_known_output(print(test_summary_3), "summary_DFOP_FOCUS_D_deSolve.txt") -}) - -context("Plotting") - -test_that("Plotting mkinfit and mmkin objects is reproducible", { - skip_on_cran() - plot_default_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]]) - plot_res_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]], show_residuals = TRUE) - plot_res_FOCUS_C_SFO_2 <- function() plot_res(fits[["SFO", "FOCUS_C"]]) - plot_sep_FOCUS_C_SFO <- function() plot_sep(fits[["SFO", "FOCUS_C"]]) - mkinparplot_FOCUS_C_SFO <- function() mkinparplot(fits[["SFO", "FOCUS_C"]]) - mkinerrplot_FOCUS_C_SFO <- function() mkinerrplot(fits[["SFO", "FOCUS_C"]]) - mmkin_FOCUS_C <- function() plot(fits[, "FOCUS_C"]) - mmkin_SFO <- function() plot(fits["SFO",]) - fit_D_obs_eigen <- suppressWarnings(mkinfit(SFO_SFO, FOCUS_2006_D, error_model = "obs", quiet = TRUE)) - fit_C_tc <- mkinfit("SFO", FOCUS_2006_C, error_model = "tc", quiet = TRUE) - - plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE) - plot_errmod_fit_C_tc <- function() plot_err(fit_C_tc) - - plot_res_sfo_sfo <- function() plot_res(f_sfo_sfo_desolve) - plot_err_sfo_sfo <- function() plot_err(f_sfo_sfo_desolve) - plot_errmod_fit_obs_1 <- function() plot_err(fit_obs_1, sep_obs = FALSE) - plot_errmod_fit_tc_1 <- function() plot_err(fit_tc_1, sep_obs = FALSE) - - skip_if(getRversion() > 4.0) - vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with defaults", plot_default_FOCUS_C_SFO) - vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with residuals like in gmkin", plot_res_FOCUS_C_SFO) - vdiffr::expect_doppelganger("plot_res for FOCUS C", plot_res_FOCUS_C_SFO_2) - vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with sep = TRUE", plot_sep_FOCUS_C_SFO) - vdiffr::expect_doppelganger("mkinparplot for FOCUS C SFO", mkinparplot_FOCUS_C_SFO) - vdiffr::expect_doppelganger("mkinerrplot for FOCUS C SFO", mkinerrplot_FOCUS_C_SFO) - vdiffr::expect_doppelganger("mmkin plot for FOCUS C", mmkin_FOCUS_C) - vdiffr::expect_doppelganger("mmkin plot for SFO (FOCUS C and D)", mmkin_SFO) - vdiffr::expect_doppelganger("plot_errmod with FOCUS C tc", plot_errmod_fit_C_tc) - skip_on_travis() # Still not working on Travis, maybe because of deSolve producing - # different results when not working with a compiled model or eigenvalues - vdiffr::expect_doppelganger("plot_errmod with FOCUS D obs eigen", plot_errmod_fit_D_obs_eigen) - vdiffr::expect_doppelganger("plot_res for FOCUS D", plot_res_sfo_sfo) - vdiffr::expect_doppelganger("plot_err for FOCUS D", plot_err_sfo_sfo) - vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_tc", plot_errmod_fit_tc_1) - vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_obs", plot_errmod_fit_obs_1) -}) - -context("AIC calculation") - -test_that("The AIC is reproducible", { - expect_equivalent(AIC(fits[["SFO", "FOCUS_C"]]), 59.3, scale = 1, tolerance = 0.1) - expect_equivalent(AIC(fits[, "FOCUS_C"]), - data.frame(df = c(3, 4, 5, 5), AIC = c(59.3, 44.7, 29.0, 39.2)), - scale = 1, tolerance = 0.1) - expect_error(AIC(fits["SFO", ]), "column object") - expect_equivalent(BIC(fits[, "FOCUS_C"]), - data.frame(df = c(3, 4, 5, 5), AIC = c(59.9, 45.5, 30.0, 40.2)), - scale = 1, tolerance = 0.1) -}) diff --git a/tests/testthat/test_summary.R b/tests/testthat/test_summary.R new file mode 100644 index 00000000..5cf6ac6b --- /dev/null +++ b/tests/testthat/test_summary.R @@ -0,0 +1,57 @@ +context("Summary") + +test_that("Summaries are reproducible", { + fit <- fits[["DFOP", "FOCUS_C"]] + test_summary <- summary(fit) + test_summary$fit_version <- "Dummy 0.0 for testing" + test_summary$fit_Rversion <- "Dummy R version for testing" + test_summary$date.fit <- "Dummy date for testing" + test_summary$date.summary <- "Dummy date for testing" + test_summary$calls <- "test 0" + test_summary$Corr <- signif(test_summary$Corr, 1) + test_summary$time <- c(elapsed = "test time 0") + # The correlation matrix is quite platform dependent + # It differs between i386 and amd64 on Windows + # and between Travis and my own Linux system + test_summary$Corr <- "Correlation matrix is platform dependent, not tested" + expect_known_output(print(test_summary), "summary_DFOP_FOCUS_C.txt") + + test_summary_2 <- summary(f_sfo_sfo_eigen) + test_summary_2$fit_version <- "Dummy 0.0 for testing" + test_summary_2$fit_Rversion <- "Dummy R version for testing" + test_summary_2$date.fit <- "Dummy date for testing" + test_summary_2$date.summary <- "Dummy date for testing" + test_summary_2$calls <- "test 0" + test_summary_2$time <- c(elapsed = "test time 0") + # The correlation matrix is quite platform dependent + # It differs between i386 and amd64 on Windows + # and between Travis and my own Linux system + # Even more so when using the Eigen method + test_summary_2$Corr <- "Correlation matrix is platform dependent, not tested" + # The residuals for this method are also platform sensitive + test_summary_2$data$residual <- "not tested" + expect_known_output(print(test_summary_2), "summary_DFOP_FOCUS_D_eigen.txt") + + test_summary_3 <- summary(f_sfo_sfo_desolve) + test_summary_3$fit_version <- "Dummy 0.0 for testing" + test_summary_3$fit_Rversion <- "Dummy R version for testing" + test_summary_3$date.fit <- "Dummy date for testing" + test_summary_3$date.summary <- "Dummy date for testing" + test_summary_3$calls <- "test 0" + test_summary_3$time <- c(elapsed = "test time 0") + # The correlation matrix is quite platform dependent + # It differs between i386 and amd64 on Windows + # and between Travis and my own Linux system + test_summary_3$Corr <- "Correlation matrix is platform dependent, not tested" + expect_known_output(print(test_summary_3), "summary_DFOP_FOCUS_D_deSolve.txt") +}) + +test_that("A fit generated with mkin 0.9.48.1 can be summarised", { + # Generated with mkin 0.9.48.1 + # SFO_SFO <- mkinmod(parent = list(type = "SFO", to = "m1"), + # m1 = list(type = "SFO"), quiet = TRUE) + # fit_old <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) + # save(fit_old, file = "~/git/mkin/inst/testdata/fit_old_FOCUS_D.rda", version = 2 ) + load(system.file("testdata/fit_old_FOCUS_D.rda", package = "mkin")) + expect_true(length(summary(fit_old)) > 0) +}) diff --git a/tests/testthat/test_twa.R b/tests/testthat/test_twa.R new file mode 100644 index 00000000..a8fc0bd9 --- /dev/null +++ b/tests/testthat/test_twa.R @@ -0,0 +1,26 @@ +context("Calculation of maximum time weighted average concentrations (TWAs)") + +test_that("Time weighted average concentrations are correct", { + skip_on_cran() + + outtimes_10 <- seq(0, 10, length.out = 10000) + + ds <- "FOCUS_C" + for (model in models) { + fit <- fits[[model, ds]] + bpar <- summary(fit)$bpar[, "Estimate"] + pred_10 <- mkinpredict(fit$mkinmod, + odeparms = bpar[2:length(bpar)], + odeini = c(parent = bpar[[1]]), + outtimes = outtimes_10) + twa_num <- mean(pred_10[, "parent"]) + names(twa_num) <- 10 + twa_ana <- max_twa_parent(fit, 10) + + # Test for absolute difference (scale = 1) + # The tolerance can be reduced if the length of outtimes is increased, + # but this needs more computing time so we stay with lenght.out = 10k + expect_equal(twa_num, twa_ana, tolerance = 0.003, scale = 1) + } +}) + |