From f0ef23a7598e5d19648ae4edc2b74e0fba27a41c Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 3 Feb 2021 16:41:31 +0100 Subject: Prepare for v1.0.0 - Improve authorship and copyright information - Prepare pkgdown config - Remove dependence on saemix as we need the development version which is not ready for CRAN - Temporarily remove saemix interface to check code coverage of the rest --- DESCRIPTION | 20 +- GNUmakefile | 2 +- NAMESPACE | 8 - NEWS.md | 15 +- R/endpoints.R | 8 +- R/plot.mixed.mmkin.R | 23 +- R/saem.R | 512 --------------------- R/summary.saem.mmkin.R | 268 ----------- _pkgdown.yml | 9 +- build.log | 3 +- check.log | 19 +- docs/404.html | 4 +- docs/articles/FOCUS_D.html | 113 +++-- docs/articles/FOCUS_D_files/figure-html/plot-1.png | Bin 81430 -> 79550 bytes docs/articles/FOCUS_L.html | 206 ++++----- .../figure-html/unnamed-chunk-10-1.png | Bin 43098 -> 42372 bytes .../figure-html/unnamed-chunk-12-1.png | Bin 84614 -> 82635 bytes .../figure-html/unnamed-chunk-13-1.png | Bin 33489 -> 32949 bytes .../figure-html/unnamed-chunk-15-1.png | Bin 59507 -> 58309 bytes .../figure-html/unnamed-chunk-4-1.png | Bin 35964 -> 35417 bytes .../figure-html/unnamed-chunk-6-1.png | Bin 36379 -> 35848 bytes .../figure-html/unnamed-chunk-8-1.png | Bin 42221 -> 41558 bytes .../figure-html/unnamed-chunk-9-1.png | Bin 41739 -> 40983 bytes docs/articles/index.html | 4 +- docs/articles/mkin.html | 18 +- .../mkin_files/figure-html/unnamed-chunk-2-1.png | Bin 92716 -> 90521 bytes docs/articles/twa.html | 6 +- docs/articles/web_only/FOCUS_Z.html | 258 +++++------ .../figure-html/FOCUS_2006_Z_fits_1-1.png | Bin 69227 -> 67049 bytes .../figure-html/FOCUS_2006_Z_fits_10-1.png | Bin 109820 -> 106768 bytes .../figure-html/FOCUS_2006_Z_fits_11-1.png | Bin 109057 -> 105777 bytes .../figure-html/FOCUS_2006_Z_fits_11a-1.png | Bin 79152 -> 76400 bytes .../figure-html/FOCUS_2006_Z_fits_11b-1.png | Bin 35736 -> 35734 bytes .../figure-html/FOCUS_2006_Z_fits_2-1.png | Bin 69227 -> 67049 bytes .../figure-html/FOCUS_2006_Z_fits_3-1.png | Bin 69006 -> 66804 bytes .../figure-html/FOCUS_2006_Z_fits_5-1.png | Bin 83658 -> 81025 bytes .../figure-html/FOCUS_2006_Z_fits_6-1.png | Bin 108929 -> 105745 bytes .../figure-html/FOCUS_2006_Z_fits_7-1.png | Bin 108277 -> 105161 bytes .../figure-html/FOCUS_2006_Z_fits_9-1.png | Bin 91586 -> 89419 bytes docs/articles/web_only/NAFTA_examples.html | 322 ++++++------- .../NAFTA_examples_files/figure-html/p10-1.png | Bin 80950 -> 79145 bytes .../NAFTA_examples_files/figure-html/p11-1.png | Bin 77737 -> 75730 bytes .../NAFTA_examples_files/figure-html/p12a-1.png | Bin 83111 -> 81215 bytes .../NAFTA_examples_files/figure-html/p12b-1.png | Bin 72088 -> 70120 bytes .../NAFTA_examples_files/figure-html/p13-1.png | Bin 79358 -> 77502 bytes .../NAFTA_examples_files/figure-html/p14-1.png | Bin 82078 -> 79961 bytes .../NAFTA_examples_files/figure-html/p15a-1.png | Bin 78259 -> 76342 bytes .../NAFTA_examples_files/figure-html/p15b-1.png | Bin 80090 -> 78251 bytes .../NAFTA_examples_files/figure-html/p16-1.png | Bin 95465 -> 93322 bytes .../NAFTA_examples_files/figure-html/p5a-1.png | Bin 84037 -> 81873 bytes .../NAFTA_examples_files/figure-html/p5b-1.png | Bin 81732 -> 79970 bytes .../NAFTA_examples_files/figure-html/p6-1.png | Bin 84055 -> 82238 bytes .../NAFTA_examples_files/figure-html/p7-1.png | Bin 103768 -> 101821 bytes .../NAFTA_examples_files/figure-html/p8-1.png | Bin 94196 -> 92053 bytes .../NAFTA_examples_files/figure-html/p9a-1.png | Bin 80132 -> 78348 bytes .../NAFTA_examples_files/figure-html/p9b-1.png | Bin 77691 -> 75683 bytes docs/articles/web_only/benchmarks.html | 118 +++-- docs/articles/web_only/compiled_models.html | 64 +-- docs/authors.html | 10 +- docs/index.html | 8 +- docs/news/index.html | 41 +- docs/pkgdown.yml | 2 +- docs/reference/AIC.mmkin.html | 4 +- docs/reference/CAKE_export.html | 46 +- docs/reference/D24_2014.html | 291 ++++++++++++ docs/reference/DFOP.solution.html | 4 +- docs/reference/Extract.mmkin.html | 76 +-- docs/reference/FOCUS_2006_DFOP_ref_A_to_B.html | 18 +- docs/reference/FOCUS_2006_FOMC_ref_A_to_F.html | 18 +- docs/reference/FOCUS_2006_HS_ref_A_to_F.html | 18 +- docs/reference/FOCUS_2006_SFO_ref_A_to_F.html | 18 +- docs/reference/FOCUS_2006_datasets.html | 28 +- docs/reference/FOMC.solution.html | 4 +- docs/reference/HS.solution.html | 4 +- docs/reference/IORE.solution.html | 4 +- docs/reference/NAFTA_SOP_2015-1.png | Bin 41379 -> 62443 bytes docs/reference/NAFTA_SOP_2015.html | 42 +- docs/reference/NAFTA_SOP_Attachment-1.png | Bin 42204 -> 64198 bytes docs/reference/NAFTA_SOP_Attachment.html | 42 +- docs/reference/Rplot001.png | Bin 1011 -> 13995 bytes docs/reference/Rplot002.png | Bin 17401 -> 13648 bytes docs/reference/Rplot003.png | Bin 15601 -> 15333 bytes docs/reference/Rplot004.png | Bin 10728 -> 10670 bytes docs/reference/Rplot005.png | Bin 19636 -> 19414 bytes docs/reference/Rplot006.png | Bin 24813 -> 24624 bytes docs/reference/Rplot007.png | Bin 25161 -> 25074 bytes docs/reference/SFO.solution.html | 4 +- docs/reference/SFORB.solution.html | 4 +- docs/reference/add_err-1.png | Bin 111359 -> 109220 bytes docs/reference/add_err-2.png | Bin 64003 -> 63336 bytes docs/reference/add_err-3.png | Bin 60019 -> 59105 bytes docs/reference/add_err.html | 12 +- docs/reference/aw.html | 34 +- docs/reference/confint.mkinfit.html | 110 ++--- docs/reference/create_deg_func.html | 30 +- docs/reference/dimethenamid_2018.html | 234 ++++++++++ docs/reference/endpoints.html | 18 +- docs/reference/experimental_data_for_UBA-1.png | Bin 107152 -> 102782 bytes docs/reference/experimental_data_for_UBA.html | 90 ++-- docs/reference/f_time_norm_focus.html | 285 ++++++++++++ docs/reference/focus_soil_moisture.html | 206 +++++++++ docs/reference/get_deg_func.html | 4 +- docs/reference/ilr.html | 4 +- docs/reference/index.html | 97 ++-- docs/reference/loftest-1.png | Bin 27354 -> 40351 bytes docs/reference/loftest-2.png | Bin 27721 -> 39817 bytes docs/reference/loftest-3.png | Bin 78754 -> 76976 bytes docs/reference/loftest-4.png | Bin 76899 -> 74678 bytes docs/reference/loftest-5.png | Bin 75806 -> 73445 bytes docs/reference/loftest.html | 93 ++-- docs/reference/logLik.mkinfit.html | 12 +- docs/reference/logistic.solution-2.png | Bin 43807 -> 42283 bytes docs/reference/logistic.solution.html | 6 +- docs/reference/lrtest.mkinfit.html | 37 +- docs/reference/max_twa_parent.html | 4 +- docs/reference/mccall81_245T-1.png | Bin 0 -> 63092 bytes docs/reference/mccall81_245T.html | 67 ++- docs/reference/mixed-1.png | Bin 0 -> 219775 bytes docs/reference/mixed.html | 279 +++++++++++ docs/reference/mkin_long_to_wide.html | 21 +- docs/reference/mkin_wide_to_long.html | 23 +- docs/reference/mkinds.html | 30 +- docs/reference/mkindsg.html | 460 ++++++++++++++++++ docs/reference/mkinerrmin.html | 33 +- docs/reference/mkinerrplot-1.png | Bin 41095 -> 41094 bytes docs/reference/mkinerrplot.html | 8 +- docs/reference/mkinfit-1.png | Bin 0 -> 66175 bytes docs/reference/mkinfit.html | 182 ++++---- docs/reference/mkinmod.html | 159 +++++-- docs/reference/mkinparplot-1.png | Bin 25704 -> 25707 bytes docs/reference/mkinparplot.html | 12 +- docs/reference/mkinplot.html | 18 +- docs/reference/mkinpredict.html | 47 +- docs/reference/mkinresplot-1.png | Bin 23814 -> 23819 bytes docs/reference/mkinresplot.html | 10 +- docs/reference/mmkin-1.png | Bin 114052 -> 111011 bytes docs/reference/mmkin-2.png | Bin 110390 -> 107610 bytes docs/reference/mmkin-3.png | Bin 97548 -> 96207 bytes docs/reference/mmkin-4.png | Bin 70005 -> 67756 bytes docs/reference/mmkin-5.png | Bin 66093 -> 65061 bytes docs/reference/mmkin.html | 63 ++- docs/reference/nafta-1.png | Bin 64837 -> 62443 bytes docs/reference/nafta.html | 6 +- docs/reference/nlme-1.png | Bin 69449 -> 70122 bytes docs/reference/nlme-2.png | Bin 0 -> 94062 bytes docs/reference/nlme.html | 39 +- docs/reference/nlme.mmkin-1.png | Bin 79142 -> 124664 bytes docs/reference/nlme.mmkin-2.png | Bin 79539 -> 169529 bytes docs/reference/nlme.mmkin.html | 343 +++++++------- docs/reference/nobs.mkinfit.html | 15 +- docs/reference/parms.html | 30 +- docs/reference/plot.mixed.mmkin-1.png | Bin 0 -> 85176 bytes docs/reference/plot.mixed.mmkin-2.png | Bin 0 -> 173876 bytes docs/reference/plot.mixed.mmkin.html | 313 +++++++++++++ docs/reference/plot.mkinfit-1.png | Bin 53731 -> 53151 bytes docs/reference/plot.mkinfit-2.png | Bin 73822 -> 73247 bytes docs/reference/plot.mkinfit-3.png | Bin 69215 -> 67810 bytes docs/reference/plot.mkinfit-4.png | Bin 73285 -> 72295 bytes docs/reference/plot.mkinfit-5.png | Bin 68646 -> 66920 bytes docs/reference/plot.mkinfit-6.png | Bin 74041 -> 73078 bytes docs/reference/plot.mkinfit-7.png | Bin 75295 -> 73909 bytes docs/reference/plot.mkinfit.html | 14 +- docs/reference/plot.mmkin-1.png | Bin 50628 -> 49402 bytes docs/reference/plot.mmkin-2.png | Bin 50911 -> 49753 bytes docs/reference/plot.mmkin-3.png | Bin 47267 -> 46090 bytes docs/reference/plot.mmkin-4.png | Bin 34181 -> 33266 bytes docs/reference/plot.mmkin-5.png | Bin 58620 -> 57777 bytes docs/reference/plot.mmkin.html | 9 +- docs/reference/plot.nafta.html | 18 +- docs/reference/print.mmkin.html | 194 ++++++++ docs/reference/reexports.html | 19 +- docs/reference/residuals.mkinfit.html | 25 +- docs/reference/schaefer07_complex_case-1.png | Bin 67741 -> 66982 bytes docs/reference/schaefer07_complex_case.html | 47 +- docs/reference/sigma_twocomp-1.png | Bin 0 -> 43910 bytes docs/reference/sigma_twocomp.html | 31 +- docs/reference/summary.mkinfit.html | 14 +- docs/reference/summary.nlme.mmkin.html | 449 ++++++++++++++++++ docs/reference/synthetic_data_for_UBA_2014-1.png | Bin 70430 -> 67756 bytes docs/reference/synthetic_data_for_UBA_2014.html | 229 ++++----- docs/reference/test_data_from_UBA_2014-1.png | Bin 59294 -> 57796 bytes docs/reference/test_data_from_UBA_2014-2.png | Bin 78194 -> 73501 bytes docs/reference/test_data_from_UBA_2014.html | 81 ++-- docs/reference/transform_odeparms.html | 171 ++++--- docs/reference/update.mkinfit-1.png | Bin 28678 -> 42360 bytes docs/reference/update.mkinfit-2.png | Bin 28541 -> 43389 bytes docs/reference/update.mkinfit.html | 29 +- docs/sitemap.xml | 35 +- man/endpoints.Rd | 8 +- man/plot.mixed.mmkin.Rd | 13 +- man/saem.Rd | 155 ------- man/summary.saem.mmkin.Rd | 100 ---- test.log | 53 +-- tests/testthat/setup_script.R | 18 - tests/testthat/summary_saem_biphasic_s.txt | 77 ---- tests/testthat/test_mixed.R | 131 ------ tests/testthat/test_plot.R | 11 - tests_slow.log | 12 +- vignettes/web_only/mkin_benchmarks.rda | Bin 1010 -> 1085 bytes 199 files changed, 4966 insertions(+), 3097 deletions(-) delete mode 100644 R/saem.R delete mode 100644 R/summary.saem.mmkin.R create mode 100644 docs/reference/D24_2014.html create mode 100644 docs/reference/dimethenamid_2018.html create mode 100644 docs/reference/f_time_norm_focus.html create mode 100644 docs/reference/focus_soil_moisture.html create mode 100644 docs/reference/mccall81_245T-1.png create mode 100644 docs/reference/mixed-1.png create mode 100644 docs/reference/mixed.html create mode 100644 docs/reference/mkindsg.html create mode 100644 docs/reference/mkinfit-1.png create mode 100644 docs/reference/nlme-2.png create mode 100644 docs/reference/plot.mixed.mmkin-1.png create mode 100644 docs/reference/plot.mixed.mmkin-2.png create mode 100644 docs/reference/plot.mixed.mmkin.html create mode 100644 docs/reference/print.mmkin.html create mode 100644 docs/reference/sigma_twocomp-1.png create mode 100644 docs/reference/summary.nlme.mmkin.html delete mode 100644 man/saem.Rd delete mode 100644 man/summary.saem.mmkin.Rd delete mode 100644 tests/testthat/summary_saem_biphasic_s.txt diff --git a/DESCRIPTION b/DESCRIPTION index fa6c1c90..8bef6055 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,14 +1,15 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 0.9.50.4 -Date: 2021-01-06 -Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), - email = "jranke@uni-bremen.de", - comment = c(ORCID = "0000-0003-4371-6538")), - person("Katrin", "Lindenberger", role = "ctb"), - person("René", "Lehmann", role = "ctb"), - person("Eurofins Regulatory AG", role = "cph")) +Version: 1.0.0 +Date: 2021-02-03 +Authors@R: c( + person("Johannes", "Ranke", role = c("aut", "cre", "cph"), + email = "jranke@uni-bremen.de", + comment = c(ORCID = "0000-0003-4371-6538")), + person("Katrin", "Lindenberger", role = "ctb", comment = "contributed to mkinresplot()"), + person("René", "Lehmann", role = "ctb", comment = "ilr() and invilr()"), + person("Eurofins Regulatory AG", role = "cph", comment = "copyright for some of the contributions of JR 2012-2014")) Description: Calculation routines based on the FOCUS Kinetics Report (2006, 2014). Includes a function for conveniently defining differential equation models, model solution based on eigenvalues if possible or using numerical @@ -20,8 +21,7 @@ Depends: R (>= 2.15.1), parallel Imports: stats, graphics, methods, deSolve, R6, inline (>= 0.3.17), numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), purrr Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, - benchmarkme, tibble, stats4, saemix (>= 3.1.9000) -Additional_repositories: https://jranke.github.io/drat + benchmarkme, tibble, stats4 License: GPL LazyLoad: yes LazyData: yes diff --git a/GNUmakefile b/GNUmakefile index aaa8cdf0..8604d3b3 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -93,7 +93,7 @@ slowtests: install vdiffr: "$(RBIN)/Rscript" -e 'vdiffr::manage_cases(filter = "plot|nafta")' -testcheck: test check +testcheck: roxygen test check README.html: README.md "$(RBIN)/Rscript" -e "rmarkdown::render('README.md', output_format = 'html_document', output_options = list(mathjax = NULL))" diff --git a/NAMESPACE b/NAMESPACE index f2497283..9776d2f3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -30,15 +30,11 @@ S3method(print,mkinmod) S3method(print,mmkin) S3method(print,nafta) S3method(print,nlme.mmkin) -S3method(print,saem.mmkin) S3method(print,summary.mkinfit) S3method(print,summary.nlme.mmkin) -S3method(print,summary.saem.mmkin) S3method(residuals,mkinfit) -S3method(saem,mmkin) S3method(summary,mkinfit) S3method(summary,nlme.mmkin) -S3method(summary,saem.mmkin) S3method(update,mkinfit) S3method(update,mmkin) S3method(update,nlme.mmkin) @@ -90,9 +86,6 @@ export(parms) export(plot_err) export(plot_res) export(plot_sep) -export(saem) -export(saemix_data) -export(saemix_model) export(sigma_twocomp) export(transform_odeparms) import(deSolve) @@ -133,6 +126,5 @@ importFrom(stats,residuals) importFrom(stats,rnorm) importFrom(stats,shapiro.test) importFrom(stats,update) -importFrom(stats,vcov) importFrom(utils,getFromNamespace) importFrom(utils,write.table) diff --git a/NEWS.md b/NEWS.md index 4e97ac01..bb3f341b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ -# mkin 0.9.50.4 (unreleased) +# mkin 1.0.0 -## General new features +## General - 'mkinmod' models gain arguments 'name' and 'dll_dir' which, in conjunction with a current version of the 'inline' package, make it possible to still use the DLL used for fast ODE solutions with 'deSolve' after saving and restoring the 'mkinmod' object. @@ -24,19 +24,10 @@ - 'mixed.mmkin' New container for mmkin objects for plotting with the 'plot.mixed.mmkin' method -- 'plot.mixed.mmkin' method used for 'nlme.mmkin' and 'saem.mmkin', both inheriting from 'mixed.mmkin' (currently virtual) +- 'plot.mixed.mmkin' method used for 'nlme.mmkin' inheriting from 'mixed.mmkin' (currently virtual) - 'plot', 'summary' and 'print' methods for 'nlme.mmkin' objects -- Add the current development version of the saemix package as a second, optional backend for mixed-effects models - -- DESCRIPTION: Additional_repositories entry pointing to my drat repository on github for a suitable saemix version - -- 'saemix_model', 'saemix_data': Helper functions to fit nonlinear mixed-effects models for mmkin row objects. - -- 'saem' generic function to fit saemix models using 'saemix_model' and 'saemix_data', with a generator 'saem.mmkin', summary and plot methods - - # mkin 0.9.50.3 (2020-10-08) - 'parms': Add a method for mmkin objects diff --git a/R/endpoints.R b/R/endpoints.R index f1f47581..b5872e68 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -10,8 +10,8 @@ #' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from #' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models #' -#' @param fit An object of class [mkinfit], [nlme.mmkin] or -#' [saem.mmkin]. Or another object that has list components +#' @param fit An object of class [mkinfit] or [nlme.mmkin] +#' or another object that has list components #' mkinmod containing an [mkinmod] degradation model, and two numeric vectors, #' bparms.optim and bparms.fixed, that contain parameter values #' for that model. @@ -20,8 +20,8 @@ #' and, if applicable, a vector of formation fractions named ff #' and, if the SFORB model was in use, a vector of eigenvalues #' of these SFORB models, equivalent to DFOP rate constants -#' @note The function is used internally by [summary.mkinfit], -#' [summary.nlme.mmkin] and [summary.saem.mmkin]. +#' @note The function is used internally by [summary.mkinfit] +#' and [summary.nlme.mmkin] #' @author Johannes Ranke #' @examples #' diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 1674d855..5a0b7412 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -2,7 +2,7 @@ utils::globalVariables("ds") #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object #' -#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin] +#' @param x An object of class [mixed.mmkin], [nlme.mmkin] #' @param i A numeric index to select datasets for which to plot the individual predictions, #' in case plots get too large #' @inheritParams plot.mkinfit @@ -39,15 +39,6 @@ utils::globalVariables("ds") #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) #' -#' f_saem <- saem(f, transformations = "saemix") -#' plot(f_saem) -#' -#' # We can overlay the two variants if we generate predictions -#' pred_nlme <- mkinpredict(dfop_sfo, -#' f_nlme$bparms.optim[-1], -#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0), -#' seq(0, 180, by = 0.2)) -#' plot(f_saem, pred_over = list(nlme = pred_nlme)) #' } #' @export plot.mixed.mmkin <- function(x, @@ -91,18 +82,6 @@ plot.mixed.mmkin <- function(x, type = ifelse(standardized, "pearson", "response")) } - if (inherits(x, "saem.mmkin")) { - if (x$transformations == "saemix") backtransform = FALSE - degparms_i <- saemix::psi(x$so) - rownames(degparms_i) <- ds_names - degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms)) - colnames(degparms_i) <- degparms_i_names - residual_type = ifelse(standardized, "standardized", "residual") - residuals <- x$data[[residual_type]] - degparms_pop <- x$so@results@fixed.effects - names(degparms_pop) <- degparms_i_names - } - degparms_fixed <- fit_1$fixed$value names(degparms_fixed) <- rownames(fit_1$fixed) degparms_all <- cbind(as.matrix(degparms_i), diff --git a/R/saem.R b/R/saem.R deleted file mode 100644 index fd2a77b4..00000000 --- a/R/saem.R +++ /dev/null @@ -1,512 +0,0 @@ -utils::globalVariables(c("predicted", "std")) - -#' Fit nonlinear mixed models with SAEM -#' -#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed -#' effects models created from [mmkin] row objects using the Stochastic Approximation -#' Expectation Maximisation algorithm (SAEM). -#' -#' An mmkin row object is essentially a list of mkinfit objects that have been -#' obtained by fitting the same model to a list of datasets using [mkinfit]. -#' -#' Starting values for the fixed effects (population mean parameters, argument -#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found -#' using [mmkin]. -#' -#' @param object An [mmkin] row object containing several fits of the same -#' [mkinmod] model to different datasets -#' @param verbose Should we print information about created objects of -#' type [saemix::SaemixModel] and [saemix::SaemixData]? -#' @param transformations Per default, all parameter transformations are done -#' in mkin. If this argument is set to 'saemix', parameter transformations -#' are done in 'saemix' for the supported cases. Currently this is only -#' supported in cases where the initial concentration of the parent is not fixed, -#' SFO or DFOP is used for the parent and there is either no metabolite or one. -#' @param degparms_start Parameter values given as a named numeric vector will -#' be used to override the starting values obtained from the 'mmkin' object. -#' @param solution_type Possibility to specify the solution type in case the -#' automatic choice is not desired -#' @param quiet Should we suppress the messages saemix prints at the beginning -#' and the end of the optimisation process? -#' @param control Passed to [saemix::saemix] -#' @param \dots Further parameters passed to [saemix::saemixModel]. -#' @return An S3 object of class 'saem.mmkin', containing the fitted -#' [saemix::SaemixObject] as a list component named 'so'. The -#' object also inherits from 'mixed.mmkin'. -#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin] -#' @examples -#' \dontrun{ -#' ds <- lapply(experimental_data_for_UBA_2019[6:10], -#' function(x) subset(x$data[c("name", "time", "value")])) -#' names(ds) <- paste("Dataset", 6:10) -#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, -#' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) -#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed) -#' -#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) -#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ]) -#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) -#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) -#' -#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use -#' # functions from saemix -#' library(saemix) -#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) -#' plot(f_saem_fomc$so, plot.type = "convergence") -#' plot(f_saem_fomc$so, plot.type = "individual.fit") -#' plot(f_saem_fomc$so, plot.type = "npde") -#' plot(f_saem_fomc$so, plot.type = "vpc") -#' -#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") -#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ]) -#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) -#' -#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), -#' A1 = mkinsub("SFO")) -#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), -#' A1 = mkinsub("SFO")) -#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), -#' A1 = mkinsub("SFO")) -#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, -#' # and compiled ODEs for FOMC that are much slower -#' f_mmkin <- mmkin(list( -#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), -#' ds, quiet = TRUE) -#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds -#' # each on this system, as we use analytical solutions written for saemix. -#' # When using the analytical solutions written for mkin this took around -#' # four minutes -#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ]) -#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ]) -#' # We can use print, plot and summary methods to check the results -#' print(f_saem_dfop_sfo) -#' plot(f_saem_dfop_sfo) -#' summary(f_saem_dfop_sfo, data = TRUE) -#' -#' # The following takes about 6 minutes -#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", -#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) -#' -#' #saemix::compare.saemix(list( -#' # f_saem_dfop_sfo$so, -#' # f_saem_dfop_sfo_deSolve$so)) -#' -#' # If the model supports it, we can also use eigenvalue based solutions, which -#' # take a similar amount of time -#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen", -#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) -#' } -#' @export -saem <- function(object, ...) UseMethod("saem") - -#' @rdname saem -#' @export -saem.mmkin <- function(object, - transformations = c("mkin", "saemix"), - degparms_start = numeric(), - solution_type = "auto", - control = list(displayProgress = FALSE, print = FALSE, - save = FALSE, save.graphs = FALSE), - verbose = FALSE, quiet = FALSE, ...) -{ - transformations <- match.arg(transformations) - m_saemix <- saemix_model(object, verbose = verbose, - degparms_start = degparms_start, solution_type = solution_type, - transformations = transformations, ...) - d_saemix <- saemix_data(object, verbose = verbose) - - fit_time <- system.time({ - utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) - }) - - transparms_optim <- f_saemix@results@fixed.effects - names(transparms_optim) <- f_saemix@results@name.fixed - - if (transformations == "mkin") { - bparms_optim <- backtransform_odeparms(transparms_optim, - object[[1]]$mkinmod, - object[[1]]$transform_rates, - object[[1]]$transform_fractions) - } else { - bparms_optim <- transparms_optim - } - - return_data <- nlme_data(object) - - return_data$predicted <- f_saemix@model@model( - psi = saemix::psi(f_saemix), - id = as.numeric(return_data$ds), - xidep = return_data[c("time", "name")]) - - return_data <- transform(return_data, - residual = predicted - value, - std = sigma_twocomp(predicted, - f_saemix@results@respar[1], f_saemix@results@respar[2])) - return_data <- transform(return_data, - standardized = residual / std) - - result <- list( - mkinmod = object[[1]]$mkinmod, - mmkin = object, - solution_type = object[[1]]$solution_type, - transformations = transformations, - transform_rates = object[[1]]$transform_rates, - transform_fractions = object[[1]]$transform_fractions, - so = f_saemix, - time = fit_time, - mean_dp_start = attr(m_saemix, "mean_dp_start"), - bparms.optim = bparms_optim, - bparms.fixed = object[[1]]$bparms.fixed, - data = return_data, - err_mod = object[[1]]$err_mod, - date.fit = date(), - saemixversion = as.character(utils::packageVersion("saemix")), - mkinversion = as.character(utils::packageVersion("mkin")), - Rversion = paste(R.version$major, R.version$minor, sep=".") - ) - - class(result) <- c("saem.mmkin", "mixed.mmkin") - return(result) -} - -#' @export -#' @rdname saem -#' @param x An saem.mmkin object to print -#' @param digits Number of digits to use for printing -print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { - cat( "Kinetic nonlinear mixed-effects model fit by SAEM" ) - cat("\nStructural model:\n") - diffs <- x$mmkin[[1]]$mkinmod$diffs - nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs) - writeLines(strwrap(nice_diffs, exdent = 11)) - cat("\nData:\n") - cat(nrow(x$data), "observations of", - length(unique(x$data$name)), "variable(s) grouped in", - length(unique(x$data$ds)), "datasets\n") - - cat("\nLikelihood computed by importance sampling\n") - print(data.frame( - AIC = AIC(x$so, type = "is"), - BIC = BIC(x$so, type = "is"), - logLik = logLik(x$so, type = "is"), - row.names = " "), digits = digits) - - cat("\nFitted parameters:\n") - conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")] - rownames(conf.int) <- x$so@results@conf.int[["name"]] - print(conf.int, digits = digits) - - invisible(x) -} - -#' @rdname saem -#' @return An [saemix::SaemixModel] object. -#' @export -saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), - degparms_start = numeric(), verbose = FALSE, ...) -{ - if (nrow(object) > 1) stop("Only row objects allowed") - - mkin_model <- object[[1]]$mkinmod - - degparms_optim <- mean_degparms(object) - if (transformations == "saemix") { - degparms_optim <- backtransform_odeparms(degparms_optim, - object[[1]]$mkinmod, - object[[1]]$transform_rates, - object[[1]]$transform_fractions) - } - degparms_fixed <- object[[1]]$bparms.fixed - - # Transformations are done in the degradation function - transform.par = rep(0, length(degparms_optim)) - - odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) - odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE) - - odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names) - odeparms_fixed <- degparms_fixed[odeparms_fixed_names] - - odeini_fixed <- degparms_fixed[odeini_fixed_parm_names] - names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names) - - model_function <- FALSE - - # Model functions with analytical solutions - # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here - # In general, we need to consider exactly how the parameters in mkinfit were specified, - # as the parameters are currently mapped by position in these solutions - sinks <- sapply(mkin_model$spec, function(x) x$sink) - if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) { - # Parent only - if (length(mkin_model$spec) == 1) { - parent_type <- mkin_model$spec[[1]]$type - if (length(odeini_fixed) == 1) { - if (parent_type == "SFO") { - stop("saemix needs at least two parameters to work on.") - } - if (parent_type == "FOMC") { - model_function <- function(psi, id, xidep) { - odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1]) - } - } - if (parent_type == "DFOP") { - model_function <- function(psi, id, xidep) { - g <- plogis(psi[id, 3]) - t <- xidep[, "time"] - odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) + - (1 - g) * exp(- exp(psi[id, 2]) * t)) - } - } - if (parent_type == "HS") { - model_function <- function(psi, id, xidep) { - tb <- exp(psi[id, 3]) - t <- xidep[, "time"] - k1 = exp(psi[id, 1]) - odeini_fixed * ifelse(t <= tb, - exp(- k1 * t), - exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb))) - } - } - } else { - if (parent_type == "SFO") { - if (transformations == "mkin") { - model_function <- function(psi, id, xidep) { - psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"]) - } - } else { - model_function <- function(psi, id, xidep) { - psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"]) - } - transform.par = c(0, 1) - } - } - if (parent_type == "FOMC") { - model_function <- function(psi, id, xidep) { - psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2]) - } - } - if (parent_type == "DFOP") { - if (transformations == "mkin") { - model_function <- function(psi, id, xidep) { - g <- plogis(psi[id, 4]) - t <- xidep[, "time"] - psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) + - (1 - g) * exp(- exp(psi[id, 3]) * t)) - } - } else { - model_function <- function(psi, id, xidep) { - g <- psi[id, 4] - t <- xidep[, "time"] - psi[id, 1] * (g * exp(- psi[id, 2] * t) + - (1 - g) * exp(- psi[id, 3] * t)) - } - transform.par = c(0, 1, 1, 3) - } - } - if (parent_type == "HS") { - model_function <- function(psi, id, xidep) { - tb <- exp(psi[id, 4]) - t <- xidep[, "time"] - k1 = exp(psi[id, 2]) - psi[id, 1] * ifelse(t <= tb, - exp(- k1 * t), - exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb))) - } - } - } - } - - # Parent with one metabolite - # Parameter names used in the model functions are as in - # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb - types <- unname(sapply(mkin_model$spec, function(x) x$type)) - if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) { - # Initial value for the metabolite (n20) must be fixed - if (names(odeini_fixed) == names(mkin_model$spec)[2]) { - n20 <- odeini_fixed - parent_name <- names(mkin_model$spec)[1] - if (identical(types, c("SFO", "SFO"))) { - if (transformations == "mkin") { - model_function <- function(psi, id, xidep) { - t <- xidep[, "time"] - n10 <- psi[id, 1] - k1 <- exp(psi[id, 2]) - k2 <- exp(psi[id, 3]) - f12 <- plogis(psi[id, 4]) - ifelse(xidep[, "name"] == parent_name, - n10 * exp(- k1 * t), - (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) + - (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1) - ) - } - } else { - model_function <- function(psi, id, xidep) { - t <- xidep[, "time"] - n10 <- psi[id, 1] - k1 <- psi[id, 2] - k2 <- psi[id, 3] - f12 <- psi[id, 4] - ifelse(xidep[, "name"] == parent_name, - n10 * exp(- k1 * t), - (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) + - (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1) - ) - } - transform.par = c(0, 1, 1, 3) - } - } - if (identical(types, c("DFOP", "SFO"))) { - if (transformations == "mkin") { - model_function <- function(psi, id, xidep) { - t <- xidep[, "time"] - n10 <- psi[id, 1] - k2 <- exp(psi[id, 2]) - f12 <- plogis(psi[id, 3]) - l1 <- exp(psi[id, 4]) - l2 <- exp(psi[id, 5]) - g <- plogis(psi[id, 6]) - ifelse(xidep[, "name"] == parent_name, - n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)), - ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) - - (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) + - ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 + - ((f12 * l1 + (f12 * g - f12) * k2) * l2 - - f12 * g * k2 * l1) * n10) * exp( - k2 * t)) / - ((l1 - k2) * l2 - k2 * l1 + k2^2) - ) - } - } else { - model_function <- function(psi, id, xidep) { - t <- xidep[, "time"] - n10 <- psi[id, 1] - k2 <- psi[id, 2] - f12 <- psi[id, 3] - l1 <- psi[id, 4] - l2 <- psi[id, 5] - g <- psi[id, 6] - ifelse(xidep[, "name"] == parent_name, - n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)), - ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) - - (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) + - ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 + - ((f12 * l1 + (f12 * g - f12) * k2) * l2 - - f12 * g * k2 * l1) * n10) * exp( - k2 * t)) / - ((l1 - k2) * l2 - k2 * l1 + k2^2) - ) - } - transform.par = c(0, 1, 3, 1, 1, 3) - } - } - } - } - } - - if (is.function(model_function) & solution_type == "auto") { - solution_type = "analytical saemix" - } else { - - if (solution_type == "auto") - solution_type <- object[[1]]$solution_type - - model_function <- function(psi, id, xidep) { - - uid <- unique(id) - - res_list <- lapply(uid, function(i) { - - transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict - names(transparms_optim) <- names(degparms_optim) - - odeini_optim <- transparms_optim[odeini_optim_parm_names] - names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names) - - odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)] - - ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names) - odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model, - transform_rates = object[[1]]$transform_rates, - transform_fractions = object[[1]]$transform_fractions) - odeparms <- c(odeparms_optim, odeparms_fixed) - - xidep_i <- subset(xidep, id == i) - - if (solution_type == "analytical") { - out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms) - } else { - - i_time <- xidep_i$time - i_name <- xidep_i$name - - out_wide <- mkinpredict(mkin_model, - odeparms = odeparms, odeini = odeini, - solution_type = solution_type, - outtimes = sort(unique(i_time)), - na_stop = FALSE - ) - - out_index <- cbind(as.character(i_time), as.character(i_name)) - out_values <- out_wide[out_index] - } - return(out_values) - }) - res <- unlist(res_list) - return(res) - } - } - - error.model <- switch(object[[1]]$err_mod, - const = "constant", - tc = "combined", - obs = "constant") - - if (object[[1]]$err_mod == "obs") { - warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model") - } - - error.init <- switch(object[[1]]$err_mod, - const = c(a = mean(sapply(object, function(x) x$errparms)), b = 1), - tc = c(a = mean(sapply(object, function(x) x$errparms[1])), - b = mean(sapply(object, function(x) x$errparms[2]))), - obs = c(a = mean(sapply(object, function(x) x$errparms)), b = 1)) - - degparms_psi0 <- degparms_optim - degparms_psi0[names(degparms_start)] <- degparms_start - psi0_matrix <- matrix(degparms_psi0, nrow = 1) - colnames(psi0_matrix) <- names(degparms_psi0) - - res <- saemix::saemixModel(model_function, - psi0 = psi0_matrix, - "Mixed model generated from mmkin object", - transform.par = transform.par, - error.model = error.model, - verbose = verbose - ) - attr(res, "mean_dp_start") <- degparms_optim - return(res) -} - -#' @rdname saem -#' @return An [saemix::SaemixData] object. -#' @export -saemix_data <- function(object, verbose = FALSE, ...) { - if (nrow(object) > 1) stop("Only row objects allowed") - ds_names <- colnames(object) - - ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")]) - names(ds_list) <- ds_names - ds_saemix_all <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") - ds_saemix <- data.frame(ds = ds_saemix_all$ds, - name = as.character(ds_saemix_all$variable), - time = ds_saemix_all$time, - value = ds_saemix_all$observed, - stringsAsFactors = FALSE) - - res <- saemix::saemixData(ds_saemix, - name.group = "ds", - name.predictors = c("time", "name"), - name.response = "value", - verbose = verbose, - ...) - return(res) -} diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R deleted file mode 100644 index e92c561c..00000000 --- a/R/summary.saem.mmkin.R +++ /dev/null @@ -1,268 +0,0 @@ -#' Summary method for class "saem.mmkin" -#' -#' Lists model equations, initial parameter values, optimised parameters -#' for fixed effects (population), random effects (deviations from the -#' population mean) and residual error model, as well as the resulting -#' endpoints such as formation fractions and DT50 values. Optionally -#' (default is FALSE), the data are listed in full. -#' -#' @param object an object of class [saem.mmkin] -#' @param x an object of class [summary.saem.mmkin] -#' @param data logical, indicating whether the full data should be included in -#' the summary. -#' @param verbose Should the summary be verbose? -#' @param distimes logical, indicating whether DT50 and DT90 values should be -#' included. -#' @param digits Number of digits to use for printing -#' @param \dots optional arguments passed to methods like \code{print}. -#' @return The summary function returns a list based on the [saemix::SaemixObject] -#' obtained in the fit, with at least the following additional components -#' \item{saemixversion, mkinversion, Rversion}{The saemix, mkin and R versions used} -#' \item{date.fit, date.summary}{The dates where the fit and the summary were -#' produced} -#' \item{diffs}{The differential equations used in the degradation model} -#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions} -#' \item{data}{The data} -#' \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals} -#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available} -#' \item{confint_errmod}{Error model parameters with confidence intervals} -#' \item{ff}{The estimated formation fractions derived from the fitted -#' model.} -#' \item{distimes}{The DT50 and DT90 values for each observed variable.} -#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} -#' The print method is called for its side effect, i.e. printing the summary. -#' @importFrom stats predict vcov -#' @author Johannes Ranke for the mkin specific parts -#' saemix authors for the parts inherited from saemix. -#' @examples -#' # Generate five datasets following DFOP-SFO kinetics -#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) -#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"), -#' m1 = mkinsub("SFO"), quiet = TRUE) -#' set.seed(1234) -#' k1_in <- rlnorm(5, log(0.1), 0.3) -#' k2_in <- rlnorm(5, log(0.02), 0.3) -#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3)) -#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3)) -#' k_m1_in <- rlnorm(5, log(0.02), 0.3) -#' -#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) { -#' mkinpredict(dfop_sfo, -#' c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1), -#' c(parent = 100, m1 = 0), -#' sampling_times) -#' } -#' -#' ds_mean_dfop_sfo <- lapply(1:5, function(i) { -#' mkinpredict(dfop_sfo, -#' c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i], -#' f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]), -#' c(parent = 100, m1 = 0), -#' sampling_times) -#' }) -#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5) -#' -#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) { -#' add_err(ds, -#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), -#' n = 1)[[1]] -#' }) -#' -#' \dontrun{ -#' # Evaluate using mmkin and saem -#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, -#' quiet = TRUE, error_model = "tc", cores = 5) -#' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo) -#' summary(f_saem_dfop_sfo, data = TRUE) -#' } -#' -#' @export -summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) { - - mod_vars <- names(object$mkinmod$diffs) - - pnames <- names(object$mean_dp_start) - np <- length(pnames) - - conf.int <- object$so@results@conf.int - rownames(conf.int) <- conf.int$name - confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")]) - colnames(confint_trans)[1] <- "est." - - # In case objects were produced by earlier versions of saem - if (is.null(object$transformations)) object$transformations <- "mkin" - - if (object$transformations == "mkin") { - bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, - object$transform_rates, object$transform_fractions) - bpnames <- names(bp) - - # Transform boundaries of CI for one parameter at a time, - # with the exception of sets of formation fractions (single fractions are OK). - f_names_skip <- character(0) - for (box in mod_vars) { # Figure out sets of fractions to skip - f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) - n_paths <- length(f_names) - if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) - } - - confint_back <- matrix(NA, nrow = length(bp), ncol = 3, - dimnames = list(bpnames, colnames(confint_trans))) - confint_back[, "est."] <- bp - - for (pname in pnames) { - if (!pname %in% f_names_skip) { - par.lower <- confint_trans[pname, "lower"] - par.upper <- confint_trans[pname, "upper"] - names(par.lower) <- names(par.upper) <- pname - bpl <- backtransform_odeparms(par.lower, object$mkinmod, - object$transform_rates, - object$transform_fractions) - bpu <- backtransform_odeparms(par.upper, object$mkinmod, - object$transform_rates, - object$transform_fractions) - confint_back[names(bpl), "lower"] <- bpl - confint_back[names(bpu), "upper"] <- bpu - } - } - } else { - confint_back <- confint_trans - } - - # Correlation of fixed effects (inspired by summary.nlme) - varFix <- vcov(object$so)[1:np, 1:np] - stdFix <- sqrt(diag(varFix)) - object$corFixed <- array( - t(varFix/stdFix)/stdFix, - dim(varFix), - list(pnames, pnames)) - - # Random effects - rnames <- paste0("SD.", pnames) - confint_ranef <- as.matrix(conf.int[rnames, c("estimate", "lower", "upper")]) - colnames(confint_ranef)[1] <- "est." - - # Error model - enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1") - confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")]) - colnames(confint_errmod)[1] <- "est." - - - object$confint_trans <- confint_trans - object$confint_ranef <- confint_ranef - object$confint_errmod <- confint_errmod - object$confint_back <- confint_back - - object$date.summary = date() - object$use_of_ff = object$mkinmod$use_of_ff - object$error_model_algorithm = object$mmkin_orig[[1]]$error_model_algorithm - err_mod = object$mmkin_orig[[1]]$err_mod - - object$diffs <- object$mkinmod$diffs - object$print_data <- data # boolean: Should we print the data? - so_pred <- object$so@results@predictions - - names(object$data)[4] <- "observed" # rename value to observed - - object$verbose <- verbose - - object$fixed <- object$mmkin_orig[[1]]$fixed - object$AIC = AIC(object$so) - object$BIC = BIC(object$so) - object$logLik = logLik(object$so, method = "is") - - ep <- endpoints(object) - if (length(ep$ff) != 0) - object$ff <- ep$ff - if (distimes) object$distimes <- ep$distimes - if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB - class(object) <- c("summary.saem.mmkin") - return(object) -} - -#' @rdname summary.saem.mmkin -#' @export -print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) { - cat("saemix version used for fitting: ", x$saemixversion, "\n") - cat("mkin version used for pre-fitting: ", x$mkinversion, "\n") - cat("R version used for fitting: ", x$Rversion, "\n") - - cat("Date of fit: ", x$date.fit, "\n") - cat("Date of summary:", x$date.summary, "\n") - - cat("\nEquations:\n") - nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) - writeLines(strwrap(nice_diffs, exdent = 11)) - - cat("\nData:\n") - cat(nrow(x$data), "observations of", - length(unique(x$data$name)), "variable(s) grouped in", - length(unique(x$data$ds)), "datasets\n") - - cat("\nModel predictions using solution type", x$solution_type, "\n") - - cat("\nFitted in", x$time[["elapsed"]], "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n") - - cat("\nVariance model: ") - cat(switch(x$err_mod, - const = "Constant variance", - obs = "Variance unique to each observed variable", - tc = "Two-component variance function"), "\n") - - cat("\nMean of starting values for individual parameters:\n") - print(x$mean_dp_start, digits = digits) - - cat("\nFixed degradation parameter values:\n") - if(length(x$fixed$value) == 0) cat("None\n") - else print(x$fixed, digits = digits) - - cat("\nResults:\n\n") - cat("Likelihood computed by importance sampling\n") - print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, - row.names = " "), digits = digits) - - cat("\nOptimised parameters:\n") - print(x$confint_trans, digits = digits) - - if (nrow(x$confint_trans) > 1) { - corr <- x$corFixed - class(corr) <- "correlation" - print(corr, title = "\nCorrelation:", ...) - } - - cat("\nRandom effects:\n") - print(x$confint_ranef, digits = digits) - - cat("\nVariance model:\n") - print(x$confint_errmod, digits = digits) - - if (x$transformations == "mkin") { - cat("\nBacktransformed parameters:\n") - print(x$confint_back, digits = digits) - } - - printSFORB <- !is.null(x$SFORB) - if(printSFORB){ - cat("\nEstimated Eigenvalues of SFORB model(s):\n") - print(x$SFORB, digits = digits,...) - } - - printff <- !is.null(x$ff) - if(printff){ - cat("\nResulting formation fractions:\n") - print(data.frame(ff = x$ff), digits = digits,...) - } - - printdistimes <- !is.null(x$distimes) - if(printdistimes){ - cat("\nEstimated disappearance times:\n") - print(x$distimes, digits = digits,...) - } - - if (x$print_data){ - cat("\nData:\n") - print(format(x$data, digits = digits, ...), row.names = FALSE) - } - - invisible(x) -} diff --git a/_pkgdown.yml b/_pkgdown.yml index 57f12fbf..7d66f5b0 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,10 +1,8 @@ url: https://pkgdown.jrwb.de/mkin -title: mkin - development: - mode: devel - version_label: info + mode: release + version_label: default template: bootswatch: spacelab @@ -42,13 +40,10 @@ reference: desc: Create and work with nonlinear mixed effects models contents: - nlme.mmkin - - saem.mmkin - plot.mixed.mmkin - summary.nlme.mmkin - - summary.saem.mmkin - nlme_function - get_deg_func - - saemix_model - mixed - title: Datasets and known results contents: diff --git a/build.log b/build.log index 6d426c06..80426597 100644 --- a/build.log +++ b/build.log @@ -6,6 +6,5 @@ * creating vignettes ... OK * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories -Removed empty directory ‘mkin/tests/testthat/_snaps’ -* building ‘mkin_0.9.50.4.tar.gz’ +* building ‘mkin_1.0.0.tar.gz’ diff --git a/check.log b/check.log index 778de34c..ba4820b3 100644 --- a/check.log +++ b/check.log @@ -5,13 +5,12 @@ * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘0.9.50.4’ +* this is package ‘mkin’ version ‘1.0.0’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke ’ * checking package namespace information ... OK -* checking package dependencies ... NOTE -Package suggested but not available for checking: ‘saemix’ +* checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for executable files ... OK @@ -46,8 +45,7 @@ Package suggested but not available for checking: ‘saemix’ * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK -* checking Rd cross-references ... NOTE -Package unavailable to check Rd xrefs: ‘saemix’ +* checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK @@ -58,10 +56,7 @@ Package unavailable to check Rd xrefs: ‘saemix’ * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... NOTE -Examples with CPU (user + system) or elapsed time > 5s - user system elapsed -D24_2014 1.808 3.425 0.986 +* checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -72,9 +67,5 @@ D24_2014 1.808 3.425 0.986 * checking for detritus in the temp directory ... OK * DONE -Status: 3 NOTEs -See - ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ -for details. - +Status: OK diff --git a/docs/404.html b/docs/404.html index 55b92696..42423c5d 100644 --- a/docs/404.html +++ b/docs/404.html @@ -71,7 +71,7 @@ mkin - 0.9.50.3 + 1.0.0 @@ -120,7 +120,7 @@
  • @@ -132,8 +135,8 @@ Compare also the code in the example section to see the degradation models." />
    • - - + +
      • @@ -174,13 +177,13 @@ Compare also the code in the example section to see the degradation models." />

      Compare also the code in the example section to see the degradation models.

      - 
      synthetic_data_for_UBA_2014
      + 
      synthetic_data_for_UBA_2014

      Format

      A list containing twelve datasets as an R6 class defined by mkinds, - each containing, among others, the following components

      + each containing, among others, the following components

      title

      The name of the dataset, e.g. SFO_lin_a

      data

      A data frame with the data in the form expected by mkinfit

      @@ -196,55 +199,59 @@ Compare also the code in the example section to see the degradation models." />

      Examples

      # \dontrun{
 # The data have been generated using the following kinetic models
-m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"),
-                           M1 = list(type = "SFO", to = "M2"),
-                           M2 = list(type = "SFO"), use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
-
-m_synth_SFO_par <- mkinmod(parent = list(type = "SFO", to = c("M1", "M2"),
-                                         sink = FALSE),
-                           M1 = list(type = "SFO"),
-                           M2 = list(type = "SFO"), use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
-m_synth_DFOP_lin <- mkinmod(parent = list(type = "DFOP", to = "M1"),
-                            M1 = list(type = "SFO", to = "M2"),
-                            M2 = list(type = "SFO"), use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
-m_synth_DFOP_par <- mkinmod(parent = list(type = "DFOP", to = c("M1", "M2"),
-                                          sink = FALSE),
-                            M1 = list(type = "SFO"),
-                            M2 = list(type = "SFO"), use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
+m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"),
+                           M1 = list(type = "SFO", to = "M2"),
+                           M2 = list(type = "SFO"), use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
+
+m_synth_SFO_par <- mkinmod(parent = list(type = "SFO", to = c("M1", "M2"),
+                                         sink = FALSE),
+                           M1 = list(type = "SFO"),
+                           M2 = list(type = "SFO"), use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
+m_synth_DFOP_lin <- mkinmod(parent = list(type = "DFOP", to = "M1"),
+                            M1 = list(type = "SFO", to = "M2"),
+                            M2 = list(type = "SFO"), use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
+m_synth_DFOP_par <- mkinmod(parent = list(type = "DFOP", to = c("M1", "M2"),
+                                          sink = FALSE),
+                            M1 = list(type = "SFO"),
+                            M2 = list(type = "SFO"), use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
 # The model predictions without intentional error were generated as follows
-sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-
-d_synth_SFO_lin <- mkinpredict(m_synth_SFO_lin,
-                               c(k_parent = 0.7, f_parent_to_M1 = 0.8,
-                                 k_M1 = 0.3, f_M1_to_M2 = 0.7,
-                                 k_M2 = 0.02),
-                               c(parent = 100, M1 = 0, M2 = 0),
-                               sampling_times)
-
-d_synth_DFOP_lin <- mkinpredict(m_synth_DFOP_lin,
-                                c(k1 = 0.2, k2 = 0.02, g = 0.5,
-                                  f_parent_to_M1 = 0.5, k_M1 = 0.3,
-                                  f_M1_to_M2 = 0.7, k_M2 = 0.02),
-                                 c(parent = 100, M1 = 0, M2 = 0),
-                                 sampling_times)
-
-d_synth_SFO_par <- mkinpredict(m_synth_SFO_par,
-                               c(k_parent = 0.2,
-                                 f_parent_to_M1 = 0.8, k_M1 = 0.01,
-                                 f_parent_to_M2 = 0.2, k_M2 = 0.02),
-                                 c(parent = 100, M1 = 0, M2 = 0),
-                                 sampling_times)
-
-d_synth_DFOP_par <- mkinpredict(m_synth_DFOP_par,
-                               c(k1 = 0.3, k2 = 0.02, g = 0.7,
-                                 f_parent_to_M1 = 0.6, k_M1 = 0.04,
-                                 f_parent_to_M2 = 0.4, k_M2 = 0.01),
-                                 c(parent = 100, M1 = 0, M2 = 0),
-                                 sampling_times)
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+
+d_synth_SFO_lin <- mkinpredict(m_synth_SFO_lin,
+                               c(k_parent = 0.7, f_parent_to_M1 = 0.8,
+                                 k_M1 = 0.3, f_M1_to_M2 = 0.7,
+                                 k_M2 = 0.02),
+                               c(parent = 100, M1 = 0, M2 = 0),
+                               sampling_times)
+
+d_synth_DFOP_lin <- mkinpredict(m_synth_DFOP_lin,
+                                c(k1 = 0.2, k2 = 0.02, g = 0.5,
+                                  f_parent_to_M1 = 0.5, k_M1 = 0.3,
+                                  f_M1_to_M2 = 0.7, k_M2 = 0.02),
+                                 c(parent = 100, M1 = 0, M2 = 0),
+                                 sampling_times)
+
+d_synth_SFO_par <- mkinpredict(m_synth_SFO_par,
+                               c(k_parent = 0.2,
+                                 f_parent_to_M1 = 0.8, k_M1 = 0.01,
+                                 f_parent_to_M2 = 0.2, k_M2 = 0.02),
+                                 c(parent = 100, M1 = 0, M2 = 0),
+                                 sampling_times)
+
+d_synth_DFOP_par <- mkinpredict(m_synth_DFOP_par,
+                               c(k1 = 0.3, k2 = 0.02, g = 0.7,
+                                 f_parent_to_M1 = 0.6, k_M1 = 0.04,
+                                 f_parent_to_M2 = 0.4, k_M2 = 0.01),
+                                 c(parent = 100, M1 = 0, M2 = 0),
+                                 sampling_times)
 
 # Construct names for datasets with errors
-d_synth_names = paste0("d_synth_", c("SFO_lin", "SFO_par",
-                                     "DFOP_lin", "DFOP_par"))
+d_synth_names = paste0("d_synth_", c("SFO_lin", "SFO_par",
+                                     "DFOP_lin", "DFOP_par"))
 
 # Original function used or adding errors. The add_err function now published
 # with this package is a slightly generalised version where the names of
@@ -265,33 +272,35 @@ Compare also the code in the example section to see the degradation models." />
 
 # The following is the simplified version of the two-component model of Rocke
 # and Lorenzato (1995)
-sdfunc_twocomp = function(value, sd_low, rsd_high) {
-  sqrt(sd_low^2 + value^2 * rsd_high^2)
-}
+sdfunc_twocomp = function(value, sd_low, rsd_high) {
+  sqrt(sd_low^2 + value^2 * rsd_high^2)
+}
 
 # Add the errors.
-for (d_synth_name in d_synth_names)
-{
-  d_synth = get(d_synth_name)
-  assign(paste0(d_synth_name, "_a"), add_err(d_synth, function(value) 3))
-  assign(paste0(d_synth_name, "_b"), add_err(d_synth, function(value) 7))
-  assign(paste0(d_synth_name, "_c"), add_err(d_synth,
-                           function(value) sdfunc_twocomp(value, 0.5, 0.07)))
+for (d_synth_name in d_synth_names)
+{
+  d_synth = get(d_synth_name)
+  assign(paste0(d_synth_name, "_a"), add_err(d_synth, function(value) 3))
+  assign(paste0(d_synth_name, "_b"), add_err(d_synth, function(value) 7))
+  assign(paste0(d_synth_name, "_c"), add_err(d_synth,
+                           function(value) sdfunc_twocomp(value, 0.5, 0.07)))
 
-}
+}
 
-d_synth_err_names = c(
-  paste(rep(d_synth_names, each = 3), letters[1:3], sep = "_")
-)
+d_synth_err_names = c(
+  paste(rep(d_synth_names, each = 3), letters[1:3], sep = "_")
+)
 
 # This is just one example of an evaluation using the kinetic model used for
 # the generation of the data
-  fit <- mkinfit(m_synth_SFO_lin, synthetic_data_for_UBA_2014[[1]]$data,
-                 quiet = TRUE)
-  plot_sep(fit)
        summary(fit)
      #> mkin version used for fitting:    0.9.50.2 
-#> R version used for fitting:       4.0.0 
-#> Date of fit:     Wed May 27 07:05:27 2020 
-#> Date of summary: Wed May 27 07:05:27 2020 
+  fit <- mkinfit(m_synth_SFO_lin, synthetic_data_for_UBA_2014[[1]]$data,
+                 quiet = TRUE)
+  plot_sep(fit)
+
        summary(fit)
+
      #> mkin version used for fitting:    1.0.0 
+#> R version used for fitting:       4.0.3 
+#> Date of fit:     Wed Feb  3 17:32:06 2021 
+#> Date of summary: Wed Feb  3 17:32:06 2021 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -300,7 +309,7 @@ Compare also the code in the example section to see the degradation models." />
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted using 817 model solutions performed in 0.627 s
+#> Fitted using 822 model solutions performed in 0.645 s
 #> 
 #> Error model: Constant variance 
 #> 
@@ -316,13 +325,13 @@ Compare also the code in the example section to see the degradation models." />
 #> f_M1_to_M2       0.5000 deparm
 #> 
 #> Starting values for the transformed parameters actually optimised:
-#>                     value lower upper
-#> parent_0       101.350000  -Inf   Inf
-#> log_k_parent    -2.302585  -Inf   Inf
-#> log_k_M1        -2.301586  -Inf   Inf
-#> log_k_M2        -2.300587  -Inf   Inf
-#> f_parent_ilr_1   0.000000  -Inf   Inf
-#> f_M1_ilr_1       0.000000  -Inf   Inf
+#>                      value lower upper
+#> parent_0        101.350000  -Inf   Inf
+#> log_k_parent     -2.302585  -Inf   Inf
+#> log_k_M1         -2.301586  -Inf   Inf
+#> log_k_M2         -2.300587  -Inf   Inf
+#> f_parent_qlogis   0.000000  -Inf   Inf
+#> f_M1_qlogis       0.000000  -Inf   Inf
 #> 
 #> Fixed parameter values:
 #>      value  type
@@ -335,32 +344,32 @@ Compare also the code in the example section to see the degradation models." />
 #>   188.7274 200.3723 -87.36368
 #> 
 #> Optimised, transformed parameters with symmetric confidence intervals:
-#>                Estimate Std. Error   Lower    Upper
-#> parent_0       102.1000    1.57000 98.8600 105.3000
-#> log_k_parent    -0.3020    0.03885 -0.3812  -0.2229
-#> log_k_M1        -1.2070    0.07123 -1.3520  -1.0620
-#> log_k_M2        -3.9010    0.06571 -4.0350  -3.7670
-#> f_parent_ilr_1   0.8492    0.16640  0.5103   1.1880
-#> f_M1_ilr_1       0.6780    0.17600  0.3196   1.0360
-#> sigma            2.2730    0.25740  1.7490   2.7970
+#>                 Estimate Std. Error   Lower    Upper
+#> parent_0        102.1000    1.57000 98.8600 105.3000
+#> log_k_parent     -0.3020    0.03885 -0.3812  -0.2229
+#> log_k_M1         -1.2070    0.07123 -1.3520  -1.0620
+#> log_k_M2         -3.9010    0.06571 -4.0350  -3.7670
+#> f_parent_qlogis   1.2010    0.23530  0.7216   1.6800
+#> f_M1_qlogis       0.9589    0.24890  0.4520   1.4660
+#> sigma             2.2730    0.25740  1.7490   2.7970
 #> 
 #> Parameter correlation:
-#>                  parent_0 log_k_parent   log_k_M1   log_k_M2 f_parent_ilr_1
-#> parent_0        1.000e+00    3.933e-01 -1.605e-01  2.819e-02     -4.624e-01
-#> log_k_parent    3.933e-01    1.000e+00 -4.082e-01  7.166e-02     -5.682e-01
-#> log_k_M1       -1.605e-01   -4.082e-01  1.000e+00 -3.929e-01      7.478e-01
-#> log_k_M2        2.819e-02    7.166e-02 -3.929e-01  1.000e+00     -2.658e-01
-#> f_parent_ilr_1 -4.624e-01   -5.682e-01  7.478e-01 -2.658e-01      1.000e+00
-#> f_M1_ilr_1      1.614e-01    4.102e-01 -8.109e-01  5.419e-01     -8.605e-01
-#> sigma          -1.384e-07   -2.581e-07  9.499e-08  1.518e-07      1.236e-07
-#>                f_M1_ilr_1      sigma
-#> parent_0        1.614e-01 -1.384e-07
-#> log_k_parent    4.102e-01 -2.581e-07
-#> log_k_M1       -8.109e-01  9.499e-08
-#> log_k_M2        5.419e-01  1.518e-07
-#> f_parent_ilr_1 -8.605e-01  1.236e-07
-#> f_M1_ilr_1      1.000e+00  8.795e-09
-#> sigma           8.795e-09  1.000e+00
+#>                   parent_0 log_k_parent   log_k_M1   log_k_M2 f_parent_qlogis
+#> parent_0         1.000e+00    3.933e-01 -1.605e-01  2.819e-02      -4.624e-01
+#> log_k_parent     3.933e-01    1.000e+00 -4.082e-01  7.166e-02      -5.682e-01
+#> log_k_M1        -1.605e-01   -4.082e-01  1.000e+00 -3.929e-01       7.478e-01
+#> log_k_M2         2.819e-02    7.166e-02 -3.929e-01  1.000e+00      -2.658e-01
+#> f_parent_qlogis -4.624e-01   -5.682e-01  7.478e-01 -2.658e-01       1.000e+00
+#> f_M1_qlogis      1.614e-01    4.102e-01 -8.109e-01  5.419e-01      -8.605e-01
+#> sigma           -7.941e-08   -9.143e-09 -1.268e-08  5.947e-08       5.657e-08
+#>                 f_M1_qlogis      sigma
+#> parent_0          1.614e-01 -7.941e-08
+#> log_k_parent      4.102e-01 -9.143e-09
+#> log_k_M1         -8.109e-01 -1.268e-08
+#> log_k_M2          5.419e-01  5.947e-08
+#> f_parent_qlogis  -8.605e-01  5.657e-08
+#> f_M1_qlogis       1.000e+00 -2.382e-10
+#> sigma            -2.382e-10  1.000e+00
 #> 
 #> Backtransformed parameters:
 #> Confidence intervals for internally transformed parameters are asymmetric.
@@ -407,8 +416,8 @@ Compare also the code in the example section to see the degradation models." />
 #>     7   parent      0.3  5.772e-01 -0.27717
 #>    14   parent      3.5  3.264e-03  3.49674
 #>    28   parent      3.2  1.045e-07  3.20000
-#>    90   parent      0.6  9.535e-10  0.60000
-#>   120   parent      3.5 -5.941e-10  3.50000
+#>    90   parent      0.6  9.532e-10  0.60000
+#>   120   parent      3.5 -5.940e-10  3.50000
 #>     1       M1     36.4  3.479e+01  1.61088
 #>     1       M1     37.4  3.479e+01  2.61088
 #>     3       M1     34.3  3.937e+01 -5.07027
@@ -418,9 +427,9 @@ Compare also the code in the example section to see the degradation models." />
 #>    14       M1      5.8  1.995e+00  3.80469
 #>    14       M1      1.2  1.995e+00 -0.79531
 #>    60       M1      0.5  2.111e-06  0.50000
-#>    90       M1      3.2 -9.676e-10  3.20000
-#>   120       M1      1.5  7.671e-10  1.50000
-#>   120       M1      0.6  7.671e-10  0.60000
+#>    90       M1      3.2 -9.671e-10  3.20000
+#>   120       M1      1.5  7.670e-10  1.50000
+#>   120       M1      0.6  7.670e-10  0.60000
 #>     1       M2      4.8  4.455e+00  0.34517
 #>     3       M2     20.9  2.153e+01 -0.62527
 #>     3       M2     19.3  2.153e+01 -2.22527
@@ -435,7 +444,7 @@ Compare also the code in the example section to see the degradation models." />
 #>    90       M2     10.6  1.013e+01  0.47130
 #>    90       M2     10.8  1.013e+01  0.67130
 #>   120       M2      9.8  5.521e+00  4.27893
-#>   120       M2      3.3  5.521e+00 -2.22107
      # }
+#>   120       M2      3.3  5.521e+00 -2.22107
      # }
      -

      Site built with pkgdown 1.5.1.

      +

      Site built with pkgdown 1.6.1.

      diff --git a/docs/reference/test_data_from_UBA_2014-1.png b/docs/reference/test_data_from_UBA_2014-1.png index 431181cf..168103ee 100644 Binary files a/docs/reference/test_data_from_UBA_2014-1.png and b/docs/reference/test_data_from_UBA_2014-1.png differ diff --git a/docs/reference/test_data_from_UBA_2014-2.png b/docs/reference/test_data_from_UBA_2014-2.png index 4a064163..68288aed 100644 Binary files a/docs/reference/test_data_from_UBA_2014-2.png and b/docs/reference/test_data_from_UBA_2014-2.png differ diff --git a/docs/reference/test_data_from_UBA_2014.html b/docs/reference/test_data_from_UBA_2014.html index 6059a4d2..eeaef9e0 100644 --- a/docs/reference/test_data_from_UBA_2014.html +++ b/docs/reference/test_data_from_UBA_2014.html @@ -73,7 +73,7 @@ mkin - 0.9.50.2 + 1.0.0 @@ -110,6 +110,9 @@
    • Example evaluation of NAFTA SOP Attachment examples
      • +
    • + Some benchmark timings +
  • @@ -118,8 +121,8 @@
    • - - + +
      • @@ -146,13 +149,13 @@ software packages (Ranke, 2014).

      - 
      test_data_from_UBA_2014
      + 
      test_data_from_UBA_2014

      Format

      A list containing three datasets as an R6 class defined by mkinds. - Each dataset has, among others, the following components

      + Each dataset has, among others, the following components

      title

      The name of the dataset, e.g. UBA_2014_WS_river

      data

      A data frame with the data in the form expected by mkinfit

      @@ -171,52 +174,62 @@ # generally considered deprecated due to the frequent occurrence of such # large parameter correlations, among other reasons (e.g. the adequacy of the # model). - m_ws <- mkinmod(parent_w = mkinsub("SFO", "parent_s"), - parent_s = mkinsub("SFO", "parent_w"))
      #> Successfully compiled differential equation model from auto-generated C code.
      f_river <- mkinfit(m_ws, test_data_from_UBA_2014[[1]]$data, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      plot_sep(f_river)
      - summary(f_river)$bpar
      #> Warning: Could not calculate correlation; no covariance matrix
      #> Estimate se_notrans t value Pr(>t) Lower Upper -#> parent_w_0 95.91998116 NA NA NA NA NA + m_ws <- mkinmod(parent_w = mkinsub("SFO", "parent_s"), + parent_s = mkinsub("SFO", "parent_w")) +
      #> Temporary DLL for differentials generated and loaded
      f_river <- mkinfit(m_ws, test_data_from_UBA_2014[[1]]$data, quiet = TRUE) +
      #> Warning: Observations with value of zero were removed from the data
      plot_sep(f_river) +
      + summary(f_river)$bpar +
      #> Warning: Could not calculate correlation; no covariance matrix
      #> Estimate se_notrans t value Pr(>t) Lower Upper +#> parent_w_0 95.91998118 NA NA NA NA NA #> k_parent_w 0.41145375 NA NA NA NA NA #> k_parent_s 0.04663944 NA NA NA NA NA #> f_parent_w_to_parent_s 0.12467894 NA NA NA NA NA #> f_parent_s_to_parent_w 0.50000000 NA NA NA NA NA -#> sigma 3.13612618 NA NA NA NA NA
      mkinerrmin(f_river)
      #> err.min n.optim df +#> sigma 3.13612618 NA NA NA NA NA
      mkinerrmin(f_river) +
      #> err.min n.optim df #> All data 0.1090929 5 6 #> parent_w 0.0817436 3 3 #> parent_s 0.1619965 2 3
      # This is the evaluation used for the validation of software packages # in the expertise from 2014 - m_soil <- mkinmod(parent = mkinsub("SFO", c("M1", "M2")), - M1 = mkinsub("SFO", "M3"), - M2 = mkinsub("SFO", "M3"), - M3 = mkinsub("SFO"), - use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      - f_soil <- mkinfit(m_soil, test_data_from_UBA_2014[[3]]$data, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      plot_sep(f_soil, lpos = c("topright", "topright", "topright", "bottomright"))
      summary(f_soil)$bpar
      #> Estimate se_notrans t value Pr(>t) Lower -#> parent_0 76.55425584 0.859186419 89.1008681 1.113862e-26 74.755958720 -#> k_parent 0.12081956 0.004601919 26.2541703 1.077361e-16 0.111561576 -#> k_M1 0.84258629 0.806165149 1.0451783 1.545282e-01 0.113778910 -#> k_M2 0.04210878 0.017083049 2.4649452 1.170195e-02 0.018013823 -#> k_M3 0.01122919 0.007245870 1.5497364 6.885076e-02 0.002909418 -#> f_parent_to_M1 0.32240193 0.240785518 1.3389590 9.819221e-02 NA -#> f_parent_to_M2 0.16099854 0.033691991 4.7785404 6.531224e-05 NA -#> f_M1_to_M3 0.27921506 0.269425582 1.0363346 1.565282e-01 0.022977955 -#> f_M2_to_M3 0.55641331 0.595121774 0.9349571 1.807710e-01 0.008002320 + m_soil <- mkinmod(parent = mkinsub("SFO", c("M1", "M2")), + M1 = mkinsub("SFO", "M3"), + M2 = mkinsub("SFO", "M3"), + M3 = mkinsub("SFO"), + use_of_ff = "max") +
      #> Temporary DLL for differentials generated and loaded
      + f_soil <- mkinfit(m_soil, test_data_from_UBA_2014[[3]]$data, quiet = TRUE) +
      #> Warning: Observations with value of zero were removed from the data
      plot_sep(f_soil, lpos = c("topright", "topright", "topright", "bottomright")) +
      summary(f_soil)$bpar +
      #> Estimate se_notrans t value Pr(>t) Lower +#> parent_0 76.55425649 0.859186399 89.1008710 1.113861e-26 74.755959406 +#> k_parent 0.12081956 0.004601918 26.2541722 1.077359e-16 0.111561575 +#> k_M1 0.84258614 0.806159820 1.0451850 1.545267e-01 0.113779670 +#> k_M2 0.04210880 0.017083035 2.4649483 1.170188e-02 0.018013857 +#> k_M3 0.01122918 0.007245855 1.5497385 6.885052e-02 0.002909431 +#> f_parent_to_M1 0.32240200 0.240783909 1.3389682 9.819073e-02 NA +#> f_parent_to_M2 0.16099855 0.033691953 4.7785463 6.531137e-05 NA +#> f_M1_to_M3 0.27921507 0.269423745 1.0363417 1.565266e-01 0.022978220 +#> f_M2_to_M3 0.55641253 0.595119954 0.9349586 1.807707e-01 0.008002509 #> sigma 1.14005399 0.149696423 7.6157731 1.727024e-07 0.826735778 #> Upper -#> parent_0 78.35255297 +#> parent_0 78.35255357 #> k_parent 0.13084582 -#> k_M1 6.23974738 -#> k_M2 0.09843271 -#> k_M3 0.04334017 +#> k_M1 6.23970352 +#> k_M2 0.09843260 +#> k_M3 0.04333992 #> f_parent_to_M1 NA #> f_parent_to_M2 NA -#> f_M1_to_M3 0.86450905 -#> f_M2_to_M3 0.99489911 -#> sigma 1.45337221
      mkinerrmin(f_soil)
      #> err.min n.optim df +#> f_M1_to_M3 0.86450768 +#> f_M2_to_M3 0.99489895 +#> sigma 1.45337221
      mkinerrmin(f_soil) +
      #> err.min n.optim df #> All data 0.09649963 9 20 #> parent 0.04721283 2 6 #> M1 0.26551208 2 5 #> M2 0.20327575 2 5 -#> M3 0.05196550 3 4
      # } +#> M3 0.05196550 3 4
      # }
      -

      Site built with pkgdown 1.5.1.

      +

      Site built with pkgdown 1.6.1.

      diff --git a/docs/reference/transform_odeparms.html b/docs/reference/transform_odeparms.html index 9b84d6bf..efaf7b46 100644 --- a/docs/reference/transform_odeparms.html +++ b/docs/reference/transform_odeparms.html @@ -77,7 +77,7 @@ the ilr transformation is used." /> mkin - 0.9.50.2 + 1.0.0 @@ -114,6 +114,9 @@ the ilr transformation is used." />
    • Example evaluation of NAFTA SOP Attachment examples
      • +
    • + Some benchmark timings +
  • @@ -122,8 +125,8 @@ the ilr transformation is used." />
    • - - + +
      • @@ -141,7 +144,7 @@ the ilr transformation is used." />

      Functions to transform and backtransform kinetic parameters for fitting

      - Source: R/transform_odeparms.R + Source: R/transform_odeparms.R
      @@ -151,22 +154,22 @@ restricted values to the full scale of real numbers. For kinetic rate constants and other parameters that can only take on positive values, a simple log transformation is used. For compositional parameters, such as the formations fractions that should always sum up to 1 and can not be negative, -the ilr transformation is used.

      +the ilr transformation is used.

      - 
      transform_odeparms(
-  parms,
-  mkinmod,
-  transform_rates = TRUE,
-  transform_fractions = TRUE
-)
+    
      transform_odeparms(
+  parms,
+  mkinmod,
+  transform_rates = TRUE,
+  transform_fractions = TRUE
+)
 
-backtransform_odeparms(
-  transparms,
-  mkinmod,
-  transform_rates = TRUE,
-  transform_fractions = TRUE
-)
      
+backtransform_odeparms(
+  transparms,
+  mkinmod,
+  transform_rates = TRUE,
+  transform_fractions = TRUE
+)

      Arguments

      @@ -178,9 +181,9 @@ equations.

      - @@ -197,10 +200,13 @@ models and the break point tb of the HS model.

      +of the estimator. The default (TRUE) is to do transformations. +The g parameter of the DFOP model is also seen as a fraction. +If a single fraction is transformed (g parameter of DFOP or only a single +target variable e.g. a single metabolite plus a pathway to sink), a +logistic transformation is used stats::qlogis(). In other cases, i.e. if +two or more formation fractions need to be transformed whose sum cannot +exceed one, the ilr transformation is used.

      @@ -211,86 +217,101 @@ fitting procedure.

      Value

      -

      A vector of transformed or backtransformed parameters with the same -names as the original parameters.

      +

      A vector of transformed or backtransformed parameters

      Details

      The transformation of sets of formation fractions is fragile, as it supposes the same ordering of the components in forward and backward transformation. -This is no problem for the internal use in mkinfit.

      -

      Functions

      +This is no problem for the internal use in mkinfit.

      +

      Author

      - -
        -

      • backtransform_odeparms: Backtransform the set of transformed parameters

        • -
      +

      Johannes Ranke

      Examples

      
-SFO_SFO <- mkinmod(
-  parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "SFO"))
      #> Successfully compiled differential equation model from auto-generated C code.
      # Fit the model to the FOCUS example dataset D using defaults
-fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      fit.s <- summary(fit)
+SFO_SFO <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = TRUE),
+  m1 = list(type = "SFO"))
+
      #> Temporary DLL for differentials generated and loaded
      # Fit the model to the FOCUS example dataset D using defaults
+fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
+
      #> Warning: Observations with value of zero were removed from the data
      fit.s <- summary(fit)
 # Transformed and backtransformed parameters
-print(fit.s$par, 3)
      #>                Estimate Std. Error   Lower   Upper
-#> parent_0         99.598     1.5702 96.4038 102.793
-#> log_k_parent     -2.316     0.0409 -2.3988  -2.233
-#> log_k_m1         -5.248     0.1332 -5.5184  -4.977
-#> f_parent_ilr_1    0.041     0.0631 -0.0875   0.169
-#> sigma             3.126     0.3585  2.3961   3.855
      print(fit.s$bpar, 3)
      #>                Estimate se_notrans t value   Pr(>t)    Lower    Upper
-#> parent_0       99.59848    1.57022   63.43 2.30e-36 96.40384 102.7931
+print(fit.s$par, 3)
+
      #>                 Estimate Std. Error  Lower  Upper
+#> parent_0         99.5985     1.5702 96.404 102.79
+#> log_k_parent     -2.3157     0.0409 -2.399  -2.23
+#> log_k_m1         -5.2475     0.1332 -5.518  -4.98
+#> f_parent_qlogis   0.0579     0.0893 -0.124   0.24
+#> sigma             3.1255     0.3585  2.396   3.85
      print(fit.s$bpar, 3)
+
      #>                Estimate se_notrans t value   Pr(>t)    Lower    Upper
+#> parent_0       99.59848    1.57022   63.43 2.30e-36 96.40383 102.7931
 #> k_parent        0.09870    0.00403   24.47 4.96e-23  0.09082   0.1073
 #> k_m1            0.00526    0.00070    7.51 6.16e-09  0.00401   0.0069
 #> f_parent_to_m1  0.51448    0.02230   23.07 3.10e-22  0.46912   0.5596
 #> sigma           3.12550    0.35852    8.72 2.24e-10  2.39609   3.8549
      
 # \dontrun{
 # Compare to the version without transforming rate parameters
-fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      #> Error in if (cost < cost.current) {    assign("cost.current", cost, inherits = TRUE)    if (!quiet)         cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),             " at call ", calls, ": ", cost.current, "\n", sep = "")}: missing value where TRUE/FALSE needed
      #> Timing stopped at: 0.002 0 0.002
      fit.2.s <- summary(fit.2)
      #> Error in summary(fit.2): object 'fit.2' not found
      print(fit.2.s$par, 3)
      #> Error in print(fit.2.s$par, 3): object 'fit.2.s' not found
      print(fit.2.s$bpar, 3)
      #> Error in print(fit.2.s$bpar, 3): object 'fit.2.s' not found
      # }
-
-initials <- fit$start$value
-names(initials) <- rownames(fit$start)
-transformed <- fit$start_transformed$value
-names(transformed) <- rownames(fit$start_transformed)
-transform_odeparms(initials, SFO_SFO)
      #>       parent_0   log_k_parent       log_k_m1 f_parent_ilr_1 
-#>     100.750000      -2.302585      -2.301586       0.000000 
      backtransform_odeparms(transformed, SFO_SFO)
      #>       parent_0       k_parent           k_m1 f_parent_to_m1 
+fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE)
+
      #> Warning: Observations with value of zero were removed from the data
      #> Error in if (cost < cost.current) {    assign("cost.current", cost, inherits = TRUE)    if (!quiet)         cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),             " at call ", calls, ": ", signif(cost.current, 6),             "\n", sep = "")}: missing value where TRUE/FALSE needed
      #> Timing stopped at: 0.003 0 0.003
      fit.2.s <- summary(fit.2)
+
      #> Error in summary(fit.2): object 'fit.2' not found
      print(fit.2.s$par, 3)
+
      #> Error in print(fit.2.s$par, 3): object 'fit.2.s' not found
      print(fit.2.s$bpar, 3)
+
      #> Error in print(fit.2.s$bpar, 3): object 'fit.2.s' not found
      # }
+
+initials <- fit$start$value
+names(initials) <- rownames(fit$start)
+transformed <- fit$start_transformed$value
+names(transformed) <- rownames(fit$start_transformed)
+transform_odeparms(initials, SFO_SFO)
+
      #>        parent_0    log_k_parent        log_k_m1 f_parent_qlogis 
+#>      100.750000       -2.302585       -2.301586        0.000000 
      backtransform_odeparms(transformed, SFO_SFO)
+
      #>       parent_0       k_parent           k_m1 f_parent_to_m1 
 #>       100.7500         0.1000         0.1001         0.5000 
      
 # \dontrun{
 # The case of formation fractions
-SFO_SFO.ff <- mkinmod(
-  parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "SFO"),
-  use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
-fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      fit.ff.s <- summary(fit.ff)
-print(fit.ff.s$par, 3)
      #>                Estimate Std. Error   Lower   Upper
-#> parent_0         99.598     1.5702 96.4038 102.793
-#> log_k_parent     -2.316     0.0409 -2.3988  -2.233
-#> log_k_m1         -5.248     0.1332 -5.5184  -4.977
-#> f_parent_ilr_1    0.041     0.0631 -0.0875   0.169
-#> sigma             3.126     0.3585  2.3961   3.855
      print(fit.ff.s$bpar, 3)
      #>                Estimate se_notrans t value   Pr(>t)    Lower    Upper
-#> parent_0       99.59848    1.57022   63.43 2.30e-36 96.40384 102.7931
+SFO_SFO.ff <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = TRUE),
+  m1 = list(type = "SFO"),
+  use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
+fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
+
      #> Warning: Observations with value of zero were removed from the data
      fit.ff.s <- summary(fit.ff)
+print(fit.ff.s$par, 3)
+
      #>                 Estimate Std. Error  Lower  Upper
+#> parent_0         99.5985     1.5702 96.404 102.79
+#> log_k_parent     -2.3157     0.0409 -2.399  -2.23
+#> log_k_m1         -5.2475     0.1332 -5.518  -4.98
+#> f_parent_qlogis   0.0579     0.0893 -0.124   0.24
+#> sigma             3.1255     0.3585  2.396   3.85
      print(fit.ff.s$bpar, 3)
+
      #>                Estimate se_notrans t value   Pr(>t)    Lower    Upper
+#> parent_0       99.59848    1.57022   63.43 2.30e-36 96.40383 102.7931
 #> k_parent        0.09870    0.00403   24.47 4.96e-23  0.09082   0.1073
 #> k_m1            0.00526    0.00070    7.51 6.16e-09  0.00401   0.0069
 #> f_parent_to_m1  0.51448    0.02230   23.07 3.10e-22  0.46912   0.5596
-#> sigma           3.12550    0.35852    8.72 2.24e-10  2.39609   3.8549
      initials <- c("f_parent_to_m1" = 0.5)
-transformed <- transform_odeparms(initials, SFO_SFO.ff)
-backtransform_odeparms(transformed, SFO_SFO.ff)
      #> f_parent_to_m1 
+#> sigma           3.12550    0.35852    8.72 2.24e-10  2.39609   3.8549
      initials <- c("f_parent_to_m1" = 0.5)
+transformed <- transform_odeparms(initials, SFO_SFO.ff)
+backtransform_odeparms(transformed, SFO_SFO.ff)
+
      #> f_parent_to_m1 
 #>            0.5 
      
 # And without sink
-SFO_SFO.ff.2 <- mkinmod(
-  parent = list(type = "SFO", to = "m1", sink = FALSE),
-  m1 = list(type = "SFO"),
-  use_of_ff = "max")
      #> Successfully compiled differential equation model from auto-generated C code.
      
-
-fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE)
      #> Warning: Observations with value of zero were removed from the data
      fit.ff.2.s <- summary(fit.ff.2)
-print(fit.ff.2.s$par, 3)
      #>              Estimate Std. Error Lower Upper
+SFO_SFO.ff.2 <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = FALSE),
+  m1 = list(type = "SFO"),
+  use_of_ff = "max")
+
      #> Temporary DLL for differentials generated and loaded
      
+
+fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE)
+
      #> Warning: Observations with value of zero were removed from the data
      fit.ff.2.s <- summary(fit.ff.2)
+print(fit.ff.2.s$par, 3)
+
      #>              Estimate Std. Error Lower Upper
 #> parent_0        84.79      3.012 78.67 90.91
 #> log_k_parent    -2.76      0.082 -2.92 -2.59
 #> log_k_m1        -4.21      0.123 -4.46 -3.96
-#> sigma            8.22      0.943  6.31 10.14
      print(fit.ff.2.s$bpar, 3)
      #>          Estimate se_notrans t value   Pr(>t)   Lower  Upper
+#> sigma            8.22      0.943  6.31 10.14
      print(fit.ff.2.s$bpar, 3)
+
      #>          Estimate se_notrans t value   Pr(>t)   Lower  Upper
 #> parent_0  84.7916    3.01203   28.15 1.92e-25 78.6704 90.913
 #> k_parent   0.0635    0.00521   12.19 2.91e-14  0.0538  0.075
 #> k_m1       0.0148    0.00182    8.13 8.81e-10  0.0115  0.019
-#> sigma      8.2229    0.94323    8.72 1.73e-10  6.3060 10.140
      # }
+#> sigma      8.2229    0.94323    8.72 1.73e-10  6.3060 10.140
      # }
      @@ -308,7 +329,7 @@ This is no problem for the internal use in mkinfit<
      -

      Site built with pkgdown 1.5.1.

      +

      Site built with pkgdown 1.6.1.

      diff --git a/docs/reference/update.mkinfit-1.png b/docs/reference/update.mkinfit-1.png index 4cbefa56..df8473c1 100644 Binary files a/docs/reference/update.mkinfit-1.png and b/docs/reference/update.mkinfit-1.png differ diff --git a/docs/reference/update.mkinfit-2.png b/docs/reference/update.mkinfit-2.png index f432f6f8..13c99b44 100644 Binary files a/docs/reference/update.mkinfit-2.png and b/docs/reference/update.mkinfit-2.png differ diff --git a/docs/reference/update.mkinfit.html b/docs/reference/update.mkinfit.html index 83b8c466..f7149d84 100644 --- a/docs/reference/update.mkinfit.html +++ b/docs/reference/update.mkinfit.html @@ -75,7 +75,7 @@ override these starting values." /> mkin - 0.9.50.2 + 1.0.0 @@ -112,6 +112,9 @@ override these starting values." />
    • Example evaluation of NAFTA SOP Attachment examples
      • +
    • + Some benchmark timings +
    • @@ -120,8 +123,8 @@ override these starting values." />
      • - - + +
        • @@ -139,7 +142,7 @@ override these starting values." />

        Update an mkinfit model with different arguments

        - Source: R/update.mkinfit.R + Source: R/update.mkinfit.R
        @@ -151,7 +154,7 @@ override these starting values.

        # S3 method for mkinfit
-update(object, ..., evaluate = TRUE)
        +update(object, ..., evaluate = TRUE)

        Arguments

      • mkinmod

      The kinetic model of class mkinmod, containing +

      The kinetic model of class mkinmod, containing the names of the model variables that are needed for grouping the -formation fractions before ilr transformation, the parameter +formation fractions before ilr transformation, the parameter names and the information if the pathway to sink is included in the model.

      Boolean specifying if formation fractions constants should be transformed in the model specification used in the fitting for better compliance with the assumption of normal distribution -of the estimator. The default (TRUE) is to do transformations. The g -parameter of the DFOP and HS models are also transformed, as they can also -be seen as compositional data. The transformation used for these -transformations is the ilr transformation.
      transparms
      @@ -175,11 +178,15 @@ remove arguments given in the original call

      Examples

      # \dontrun{
-fit <- mkinfit("SFO", subset(FOCUS_2006_D, value != 0), quiet = TRUE)
-parms(fit)
      #>      parent_0 k_parent_sink         sigma 
-#>   99.44423886    0.09793574    3.39632469 
      plot_err(fit)
      fit_2 <- update(fit, error_model = "tc")
-parms(fit_2)
      #>      parent_0 k_parent_sink     sigma_low      rsd_high 
-#>  1.008549e+02  1.005665e-01  3.752222e-03  6.763434e-02 
      plot_err(fit_2)
      # }
+fit <- mkinfit("SFO", subset(FOCUS_2006_D, value != 0), quiet = TRUE)
+parms(fit)
+
      #>    parent_0    k_parent       sigma 
+#> 99.44423885  0.09793574  3.39632469 
      plot_err(fit)
+
      fit_2 <- update(fit, error_model = "tc")
+parms(fit_2)
+
      #>     parent_0     k_parent    sigma_low     rsd_high 
+#> 1.008549e+02 1.005665e-01 3.752222e-03 6.763434e-02 
      plot_err(fit_2)
+
      # }
      -

      Site built with pkgdown 1.5.1.

      +

      Site built with pkgdown 1.6.1.

      diff --git a/docs/sitemap.xml b/docs/sitemap.xml index 81368436..7e3d8a53 100644 --- a/docs/sitemap.xml +++ b/docs/sitemap.xml @@ -9,6 +9,9 @@ https://pkgdown.jrwb.de/mkin/reference/CAKE_export.html + + https://pkgdown.jrwb.de/mkin/reference/D24_2014.html + https://pkgdown.jrwb.de/mkin/reference/DFOP.solution.html @@ -63,12 +66,21 @@ https://pkgdown.jrwb.de/mkin/reference/create_deg_func.html + + https://pkgdown.jrwb.de/mkin/reference/dimethenamid_2018.html + https://pkgdown.jrwb.de/mkin/reference/endpoints.html https://pkgdown.jrwb.de/mkin/reference/experimental_data_for_UBA.html + + https://pkgdown.jrwb.de/mkin/reference/f_time_norm_focus.html + + + https://pkgdown.jrwb.de/mkin/reference/focus_soil_moisture.html + https://pkgdown.jrwb.de/mkin/reference/get_deg_func.html @@ -93,6 +105,9 @@ https://pkgdown.jrwb.de/mkin/reference/mccall81_245T.html + + https://pkgdown.jrwb.de/mkin/reference/mixed.html + https://pkgdown.jrwb.de/mkin/reference/mkin_long_to_wide.html @@ -102,6 +117,9 @@ https://pkgdown.jrwb.de/mkin/reference/mkinds.html + + https://pkgdown.jrwb.de/mkin/reference/mkindsg.html + https://pkgdown.jrwb.de/mkin/reference/mkinerrmin.html @@ -126,9 +144,6 @@ https://pkgdown.jrwb.de/mkin/reference/mkinresplot.html - - https://pkgdown.jrwb.de/mkin/reference/mkinsub.html - https://pkgdown.jrwb.de/mkin/reference/mmkin.html @@ -147,6 +162,9 @@ https://pkgdown.jrwb.de/mkin/reference/parms.html + + https://pkgdown.jrwb.de/mkin/reference/plot.mixed.mmkin.html + https://pkgdown.jrwb.de/mkin/reference/plot.mkinfit.html @@ -157,13 +175,7 @@ https://pkgdown.jrwb.de/mkin/reference/plot.nafta.html - https://pkgdown.jrwb.de/mkin/reference/plot.nlme.mmkin.html - - - https://pkgdown.jrwb.de/mkin/reference/print.mkinds.html - - - https://pkgdown.jrwb.de/mkin/reference/print.mkinmod.html + https://pkgdown.jrwb.de/mkin/reference/print.mmkin.html https://pkgdown.jrwb.de/mkin/reference/reexports.html @@ -180,6 +192,9 @@ https://pkgdown.jrwb.de/mkin/reference/summary.mkinfit.html + + https://pkgdown.jrwb.de/mkin/reference/summary.nlme.mmkin.html + https://pkgdown.jrwb.de/mkin/reference/synthetic_data_for_UBA_2014.html diff --git a/man/endpoints.Rd b/man/endpoints.Rd index 72487717..0b225e62 100644 --- a/man/endpoints.Rd +++ b/man/endpoints.Rd @@ -8,8 +8,8 @@ with mkinfit} endpoints(fit) } \arguments{ -\item{fit}{An object of class \link{mkinfit}, \link{nlme.mmkin} or -\link{saem.mmkin}. Or another object that has list components +\item{fit}{An object of class \link{mkinfit} or \link{nlme.mmkin} +or another object that has list components mkinmod containing an \link{mkinmod} degradation model, and two numeric vectors, bparms.optim and bparms.fixed, that contain parameter values for that model.} @@ -32,8 +32,8 @@ Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models } \note{ -The function is used internally by \link{summary.mkinfit}, -\link{summary.nlme.mmkin} and \link{summary.saem.mmkin}. +The function is used internally by \link{summary.mkinfit} +and \link{summary.nlme.mmkin} } \examples{ diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd index c7b2344f..87a82286 100644 --- a/man/plot.mixed.mmkin.Rd +++ b/man/plot.mixed.mmkin.Rd @@ -17,7 +17,7 @@ maxabs = "auto", ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)), nrow.legend = ceiling((length(i) + 1)/ncol.legend), - rel.height.legend = 0.03 + 0.08 * nrow.legend, + rel.height.legend = 0.02 + 0.07 * nrow.legend, rel.height.bottom = 1.1, pch_ds = 1:length(i), col_ds = pch_ds + 1, @@ -27,7 +27,7 @@ ) } \arguments{ -\item{x}{An object of class \link{mixed.mmkin}, \link{saem.mmkin} or \link{nlme.mmkin}} +\item{x}{An object of class \link{mixed.mmkin}, \link{nlme.mmkin}} \item{i}{A numeric index to select datasets for which to plot the individual predictions, in case plots get too large} @@ -94,15 +94,6 @@ plot(f[, 3:4], standardized = TRUE) f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) plot(f_nlme) -f_saem <- saem(f, transformations = "saemix") -plot(f_saem) - -# We can overlay the two variants if we generate predictions -pred_nlme <- mkinpredict(dfop_sfo, - f_nlme$bparms.optim[-1], - c(parent = f_nlme$bparms.optim[[1]], A1 = 0), - seq(0, 180, by = 0.2)) -plot(f_saem, pred_over = list(nlme = pred_nlme)) } } \author{ diff --git a/man/saem.Rd b/man/saem.Rd deleted file mode 100644 index d5a8f17e..00000000 --- a/man/saem.Rd +++ /dev/null @@ -1,155 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/saem.R -\name{saem} -\alias{saem} -\alias{saem.mmkin} -\alias{print.saem.mmkin} -\alias{saemix_model} -\alias{saemix_data} -\title{Fit nonlinear mixed models with SAEM} -\usage{ -saem(object, ...) - -\method{saem}{mmkin}( - object, - transformations = c("mkin", "saemix"), - degparms_start = numeric(), - solution_type = "auto", - control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = - FALSE), - verbose = FALSE, - quiet = FALSE, - ... -) - -\method{print}{saem.mmkin}(x, digits = max(3, getOption("digits") - 3), ...) - -saemix_model( - object, - solution_type = "auto", - transformations = c("mkin", "saemix"), - degparms_start = numeric(), - verbose = FALSE, - ... -) - -saemix_data(object, verbose = FALSE, ...) -} -\arguments{ -\item{object}{An \link{mmkin} row object containing several fits of the same -\link{mkinmod} model to different datasets} - -\item{\dots}{Further parameters passed to \link[saemix:saemixModel]{saemix::saemixModel}.} - -\item{transformations}{Per default, all parameter transformations are done -in mkin. If this argument is set to 'saemix', parameter transformations -are done in 'saemix' for the supported cases. Currently this is only -supported in cases where the initial concentration of the parent is not fixed, -SFO or DFOP is used for the parent and there is either no metabolite or one.} - -\item{degparms_start}{Parameter values given as a named numeric vector will -be used to override the starting values obtained from the 'mmkin' object.} - -\item{solution_type}{Possibility to specify the solution type in case the -automatic choice is not desired} - -\item{control}{Passed to \link[saemix:saemix]{saemix::saemix}} - -\item{verbose}{Should we print information about created objects of -type \link[saemix:SaemixModel-class]{saemix::SaemixModel} and \link[saemix:SaemixData-class]{saemix::SaemixData}?} - -\item{quiet}{Should we suppress the messages saemix prints at the beginning -and the end of the optimisation process?} - -\item{x}{An saem.mmkin object to print} - -\item{digits}{Number of digits to use for printing} -} -\value{ -An S3 object of class 'saem.mmkin', containing the fitted -\link[saemix:SaemixObject-class]{saemix::SaemixObject} as a list component named 'so'. The -object also inherits from 'mixed.mmkin'. - -An \link[saemix:SaemixModel-class]{saemix::SaemixModel} object. - -An \link[saemix:SaemixData-class]{saemix::SaemixData} object. -} -\description{ -This function uses \code{\link[saemix:saemix]{saemix::saemix()}} as a backend for fitting nonlinear mixed -effects models created from \link{mmkin} row objects using the Stochastic Approximation -Expectation Maximisation algorithm (SAEM). -} -\details{ -An mmkin row object is essentially a list of mkinfit objects that have been -obtained by fitting the same model to a list of datasets using \link{mkinfit}. - -Starting values for the fixed effects (population mean parameters, argument -psi0 of \code{\link[saemix:saemixModel]{saemix::saemixModel()}} are the mean values of the parameters found -using \link{mmkin}. -} -\examples{ -\dontrun{ -ds <- lapply(experimental_data_for_UBA_2019[6:10], - function(x) subset(x$data[c("name", "time", "value")])) -names(ds) <- paste("Dataset", 6:10) -f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, - state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) -f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed) - -f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) -f_saem_sfo <- saem(f_mmkin_parent["SFO", ]) -f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) -f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) - -# The returned saem.mmkin object contains an SaemixObject, therefore we can use -# functions from saemix -library(saemix) -compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) -plot(f_saem_fomc$so, plot.type = "convergence") -plot(f_saem_fomc$so, plot.type = "individual.fit") -plot(f_saem_fomc$so, plot.type = "npde") -plot(f_saem_fomc$so, plot.type = "vpc") - -f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") -f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ]) -compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) - -sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), - A1 = mkinsub("SFO")) -fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), - A1 = mkinsub("SFO")) -dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), - A1 = mkinsub("SFO")) -# The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, -# and compiled ODEs for FOMC that are much slower -f_mmkin <- mmkin(list( - "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), - ds, quiet = TRUE) -# saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds -# each on this system, as we use analytical solutions written for saemix. -# When using the analytical solutions written for mkin this took around -# four minutes -f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ]) -f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ]) -# We can use print, plot and summary methods to check the results -print(f_saem_dfop_sfo) -plot(f_saem_dfop_sfo) -summary(f_saem_dfop_sfo, data = TRUE) - -# The following takes about 6 minutes -#f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", -# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) - -#saemix::compare.saemix(list( -# f_saem_dfop_sfo$so, -# f_saem_dfop_sfo_deSolve$so)) - -# If the model supports it, we can also use eigenvalue based solutions, which -# take a similar amount of time -#f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen", -# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) -} -} -\seealso{ -\link{summary.saem.mmkin} \link{plot.mixed.mmkin} -} diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd deleted file mode 100644 index 67cb3cbb..00000000 --- a/man/summary.saem.mmkin.Rd +++ /dev/null @@ -1,100 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/summary.saem.mmkin.R -\name{summary.saem.mmkin} -\alias{summary.saem.mmkin} -\alias{print.summary.saem.mmkin} -\title{Summary method for class "saem.mmkin"} -\usage{ -\method{summary}{saem.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) - -\method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) -} -\arguments{ -\item{object}{an object of class \link{saem.mmkin}} - -\item{data}{logical, indicating whether the full data should be included in -the summary.} - -\item{verbose}{Should the summary be verbose?} - -\item{distimes}{logical, indicating whether DT50 and DT90 values should be -included.} - -\item{\dots}{optional arguments passed to methods like \code{print}.} - -\item{x}{an object of class \link{summary.saem.mmkin}} - -\item{digits}{Number of digits to use for printing} -} -\value{ -The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject} -obtained in the fit, with at least the following additional components -\item{saemixversion, mkinversion, Rversion}{The saemix, mkin and R versions used} -\item{date.fit, date.summary}{The dates where the fit and the summary were -produced} -\item{diffs}{The differential equations used in the degradation model} -\item{use_of_ff}{Was maximum or minimum use made of formation fractions} -\item{data}{The data} -\item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals} -\item{confint_back}{Backtransformed parameters, with confidence intervals if available} -\item{confint_errmod}{Error model parameters with confidence intervals} -\item{ff}{The estimated formation fractions derived from the fitted -model.} -\item{distimes}{The DT50 and DT90 values for each observed variable.} -\item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} -The print method is called for its side effect, i.e. printing the summary. -} -\description{ -Lists model equations, initial parameter values, optimised parameters -for fixed effects (population), random effects (deviations from the -population mean) and residual error model, as well as the resulting -endpoints such as formation fractions and DT50 values. Optionally -(default is FALSE), the data are listed in full. -} -\examples{ -# Generate five datasets following DFOP-SFO kinetics -sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) -dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"), - m1 = mkinsub("SFO"), quiet = TRUE) -set.seed(1234) -k1_in <- rlnorm(5, log(0.1), 0.3) -k2_in <- rlnorm(5, log(0.02), 0.3) -g_in <- plogis(rnorm(5, qlogis(0.5), 0.3)) -f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3)) -k_m1_in <- rlnorm(5, log(0.02), 0.3) - -pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) { - mkinpredict(dfop_sfo, - c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1), - c(parent = 100, m1 = 0), - sampling_times) -} - -ds_mean_dfop_sfo <- lapply(1:5, function(i) { - mkinpredict(dfop_sfo, - c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i], - f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]), - c(parent = 100, m1 = 0), - sampling_times) -}) -names(ds_mean_dfop_sfo) <- paste("ds", 1:5) - -ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) { - add_err(ds, - sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), - n = 1)[[1]] -}) - -\dontrun{ -# Evaluate using mmkin and saem -f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, - quiet = TRUE, error_model = "tc", cores = 5) -f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo) -summary(f_saem_dfop_sfo, data = TRUE) -} - -} -\author{ -Johannes Ranke for the mkin specific parts -saemix authors for the parts inherited from saemix. -} diff --git a/test.log b/test.log index ead5a1c7..1e64f466 100644 --- a/test.log +++ b/test.log @@ -4,47 +4,34 @@ Testing mkin ✔ | OK F W S | Context ✔ | 5 | AIC calculation ✔ | 2 | Export dataset for reading into CAKE -✔ | 14 | 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) [1.0 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s] -✔ | 7 | Fitting the SFORB model [4.0 s] -✔ | 5 | Analytical solutions for coupled models [3.6 s] +✔ | 7 | Fitting the SFORB model [3.5 s] +✔ | 5 | Analytical solutions for coupled models [3.3 s] ✔ | 5 | Calculation of Akaike weights -✔ | 14 | Confidence intervals and p-values [1.4 s] -✔ | 14 | Error model fitting [6.8 s] +✔ | 14 | Confidence intervals and p-values [1.2 s] +✔ | 14 | Error model fitting [4.4 s] ✔ | 5 | Time step normalisation -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.5 s] -✔ | 1 | Fitting the logistic model [0.3 s] -✔ | 5 3 | Nonlinear mixed-effects models [0.1 s] -──────────────────────────────────────────────────────────────────────────────── -Skip (test_mixed.R:4:3): Parent fits using saemix are correctly implemented -Reason: !saemix_available is TRUE - -Skip (test_mixed.R:94:3): Print methods work -Reason: !saemix_available is TRUE - -Skip (test_mixed.R:118:3): saem results are reproducible for biphasic fits -Reason: !saemix_available is TRUE -──────────────────────────────────────────────────────────────────────────────── +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] +✔ | 1 | Fitting the logistic model [0.2 s] +✔ | 5 | Nonlinear mixed-effects models [0.1 s] ✔ | 2 | Test dataset classes mkinds and mkindsg -✔ | 1 | mkinfit features [0.4 s] -✔ | 10 | Special cases of mkinfit calls [0.4 s] +✔ | 1 | mkinfit features [0.3 s] +✔ | 10 | Special cases of mkinfit calls [0.3 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 [2.0 s] -✔ | 9 | Nonlinear mixed-effects models [9.3 s] -✔ | 14 | Plotting [2.0 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.6 s] +✔ | 9 | Nonlinear mixed-effects models [7.7 s] +✔ | 14 | Plotting [1.7 s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.8 s] -✔ | 4 | Summary [0.2 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] +✔ | 4 | Summary [0.1 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.5 s] -✔ | 9 | Hypothesis tests [7.9 s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.7 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] +✔ | 9 | Hypothesis tests [7.2 s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 48.3 s - -── Skipped tests ────────────────────────────────────────────────────────────── -● !saemix_available is TRUE (3) +Duration: 40.4 s -[ FAIL 0 | WARN 0 | SKIP 3 | PASS 176 ] +[ FAIL 0 | WARN 0 | SKIP 0 | PASS 176 ] diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 2071e05c..547b2d6c 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -178,10 +178,6 @@ ds_biphasic <- lapply(ds_biphasic_mean, function(ds) { }) # Mixed model fits -saemix_available <- FALSE -if (requireNamespace("saemix", quietly = TRUE)) { - if(packageVersion("saemix") > "3.1.9000") saemix_available <- TRUE -} mmkin_sfo_1 <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = n_cores) mmkin_dfop_1 <- mmkin("DFOP", ds_dfop, quiet = TRUE, cores = n_cores) mmkin_biphasic <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, quiet = TRUE, cores = n_cores) @@ -190,16 +186,6 @@ mmkin_biphasic_mixed <- mixed(mmkin_biphasic) dfop_nlme_1 <- nlme(mmkin_dfop_1) nlme_biphasic <- nlme(mmkin_biphasic) -if (saemix_available) { - sfo_saem_1 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "saemix") - - dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin") - dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix") - - saem_biphasic_m <- saem(mmkin_biphasic, transformations = "mkin", quiet = TRUE) - saem_biphasic_s <- saem(mmkin_biphasic, transformations = "saemix", quiet = TRUE) -} - ds_uba <- lapply(experimental_data_for_UBA_2019[6:10], function(x) subset(x$data[c("name", "time", "value")])) names(ds_uba) <- paste("Dataset", 6:10) @@ -211,7 +197,3 @@ f_uba_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo_uba, "DFOP-SFO" = dfop_sfo_uba), ds_uba, quiet = TRUE, cores = n_cores) f_uba_dfop_sfo_mixed <- mixed(f_uba_mmkin[2, ]) -if (saemix_available) { - f_uba_sfo_sfo_saem <- saem(f_uba_mmkin["SFO-SFO", ], quiet = TRUE, transformations = "saemix") - f_uba_dfop_sfo_saem <- saem(f_uba_mmkin["DFOP-SFO", ], quiet = TRUE, transformations = "saemix") -} diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt deleted file mode 100644 index b5777ce1..00000000 --- a/tests/testthat/summary_saem_biphasic_s.txt +++ /dev/null @@ -1,77 +0,0 @@ -saemix version used for fitting: Dummy 0.0 for testing -mkin version used for pre-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 -d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) - * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * - exp(-k2 * time))) * parent - k_m1 * m1 - -Data: -509 observations of 2 variable(s) grouped in 15 datasets - -Model predictions using solution type analytical - -Fitted in test time 0 s using 300, 100 iterations - -Variance model: Constant variance - -Mean of starting values for individual parameters: - parent_0 k_m1 f_parent_to_m1 k1 k2 - 1.0e+02 1.8e-03 4.8e-01 4.5e-02 1.9e-02 - g - 5.1e-01 - -Fixed degradation parameter values: -None - -Results: - -Likelihood computed by importance sampling - AIC BIC logLik - 2692 2702 -1333 - -Optimised parameters: - est. lower upper -parent_0 1.0e+02 99.3030 1.0e+02 -k_m1 2.2e-03 0.0016 2.8e-03 -f_parent_to_m1 4.9e-01 0.4477 5.3e-01 -k1 4.0e-02 0.0239 5.6e-02 -k2 1.5e-02 0.0085 2.2e-02 -g 5.0e-01 0.3805 6.2e-01 - -Correlation: - prnt_0 k_m1 f_p__1 k1 k2 -k_m1 -0.152 -f_parent_to_m1 -0.138 0.396 -k1 0.085 -0.036 -0.026 -k2 0.005 0.067 0.018 -0.043 -g -0.075 -0.113 -0.029 -0.248 -0.263 - -Random effects: - est. lower upper -SD.parent_0 1.23 0.39 2.07 -SD.k_m1 0.14 -0.48 0.75 -SD.f_parent_to_m1 0.28 0.17 0.39 -SD.k1 0.71 0.42 1.00 -SD.k2 0.78 0.47 1.09 -SD.g 0.27 -0.33 0.86 - -Variance model: - est. lower upper -a.1 2.7 2.6 2.9 - -Resulting formation fractions: - ff -parent_m1 0.49 -parent_sink 0.51 - -Estimated disappearance times: - DT50 DT90 DT50back DT50_k1 DT50_k2 -parent 27 110 33 17 46 -m1 318 1055 NA NA NA diff --git a/tests/testthat/test_mixed.R b/tests/testthat/test_mixed.R index e9af10e6..ca0072ef 100644 --- a/tests/testthat/test_mixed.R +++ b/tests/testthat/test_mixed.R @@ -1,98 +1,9 @@ context("Nonlinear mixed-effects models") -test_that("Parent fits using saemix are correctly implemented", { - skip_if(!saemix_available) - - expect_error(saem(fits), "Only row objects") - # Some fits were done in the setup script - mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100)) - expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*") - - sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin") - sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters") - s_sfo_s1 <- summary(sfo_saem_1) - s_sfo_s2 <- summary(sfo_saem_2) - - sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge") - s_sfo_n <- summary(sfo_nlme_1) - - # Compare with input - expect_equal(round(s_sfo_s2$confint_ranef["SD.log_k_parent", "est."], 1), 0.3) - # k_parent is a bit different from input 0.03 here - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035) - expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035) - - # But the result is pretty unanimous between methods - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), - round(s_sfo_s2$confint_back["k_parent", "est."], 3)) - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), - round(s_sfo_n$confint_back["k_parent", "est."], 3)) - - mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores) - fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE) - ci_fomc_s1 <- summary(fomc_saem_1)$confint_back - - fomc_pop <- as.numeric(fomc_pop) - expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop)) - expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop)) - - mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent") - fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin") - ci_fomc_s2 <- summary(fomc_saem_2)$confint_back - - expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) - expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) - - s_dfop_s1 <- summary(dfop_saemix_1) - s_dfop_s2 <- summary(dfop_saemix_2) - s_dfop_n <- summary(dfop_nlme_1) - - dfop_pop <- as.numeric(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)) - - dfop_mmkin_means_trans <- apply(parms(mmkin_dfop_1, transformed = TRUE), 1, mean) - dfop_mmkin_means <- backtransform_odeparms(dfop_mmkin_means_trans, mmkin_dfop_1$mkinmod) - - # We get < 22% deviations by averaging the transformed parameters - rel_diff_mmkin <- (dfop_mmkin_means - dfop_pop) / dfop_pop - expect_true(all(rel_diff_mmkin < 0.22)) - - # We get < 50% 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.5)) - - # We get < 12% deviations with transformations made in saemix - rel_diff_2 <- (s_dfop_s2$confint_back[, "est."] - dfop_pop) / dfop_pop - expect_true(all(rel_diff_2 < 0.12)) - - mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) - hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) - 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_2 <- saem(mmkin_hs_2, quiet = TRUE) - ci_hs_s2 <- summary(hs_saem_2)$confint_back - - #expect_true(all(ci_hs_s2[, "lower"] < hs_pop[2:4])) # k1 again overestimated - expect_true(all(ci_hs_s2[, "upper"] > hs_pop[2:4])) - - # HS would likely benefit from implemenation of transformations = "saemix" -}) - test_that("Print methods work", { expect_known_output(print(fits, digits = 2), "print_mmkin_parent.txt") expect_known_output(print(mmkin_biphasic_mixed, digits = 2), "print_mmkin_biphasic_mixed.txt") expect_known_output(print(nlme_biphasic, digits = 1), "print_nlme_biphasic.txt") - - skip_if(!saemix_available) - expect_known_output(print(sfo_saem_1, digits = 1), "print_sfo_saem_1.txt") }) test_that("nlme results are reproducible to some degree", { @@ -112,45 +23,3 @@ test_that("nlme results are reproducible to some degree", { # expect_true(all(ci_dfop_sfo_n[, "lower"] < dfop_sfo_pop)) # k2 is overestimated expect_true(all(ci_dfop_sfo_n[, "upper"] > dfop_sfo_pop)) }) - -test_that("saem results are reproducible for biphasic fits", { - - skip_if(!saemix_available) - test_summary <- summary(saem_biphasic_s) - test_summary$saemixversion <- "Dummy 0.0 for testing" - test_summary$mkinversion <- "Dummy 0.0 for testing" - test_summary$Rversion <- "Dummy R version for testing" - test_summary$date.fit <- "Dummy date for testing" - test_summary$date.summary <- "Dummy date for testing" - test_summary$time <- c(elapsed = "test time 0") - - expect_known_output(print(test_summary, digits = 2), "summary_saem_biphasic_s.txt") - - dfop_sfo_pop <- as.numeric(dfop_sfo_pop) - no_k1 <- c(1, 2, 3, 5, 6) - no_k2 <- c(1, 2, 3, 4, 6) - - ci_dfop_sfo_s_s <- summary(saem_biphasic_s)$confint_back - expect_true(all(ci_dfop_sfo_s_s[, "lower"] < dfop_sfo_pop)) - expect_true(all(ci_dfop_sfo_s_s[, "upper"] > dfop_sfo_pop)) - - # k1 and k2 are not fitted well - ci_dfop_sfo_s_m <- summary(saem_biphasic_m)$confint_back - expect_true(all(ci_dfop_sfo_s_m[no_k2, "lower"] < dfop_sfo_pop[no_k2])) - expect_true(all(ci_dfop_sfo_s_m[no_k1, "upper"] > dfop_sfo_pop[no_k1])) - - # I tried to only do few iterations in routine tests as this is so slow - # but then deSolve fails at some point (presumably at the switch between - # the two types of iterations) - #saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve", - # control = list(nbiter.saemix = c(10, 5), nbiter.burn = 5), quiet = TRUE) - - skip("Fitting with saemix takes around 10 minutes when using deSolve") - saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve", quiet = TRUE) - - # As with the analytical solution, k1 and k2 are not fitted well - ci_dfop_sfo_s_d <- summary(saem_biphasic_2)$confint_back - expect_true(all(ci_dfop_sfo_s_d[no_k2, "lower"] < dfop_sfo_pop[no_k2])) - expect_true(all(ci_dfop_sfo_s_d[no_k1, "upper"] > dfop_sfo_pop[no_k1])) -}) - diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R index 1c95d069..0bf3ee66 100644 --- a/tests/testthat/test_plot.R +++ b/tests/testthat/test_plot.R @@ -35,11 +35,6 @@ test_that("Plotting mkinfit, mmkin and mixed model objects is reproducible", { plot_biphasic_mmkin <- function() plot(f_uba_dfop_sfo_mixed) vdiffr::expect_doppelganger("mixed model fit for mmkin object", plot_biphasic_mmkin) - if (saemix_available) { - plot_biphasic_saem_s <- function() plot(f_uba_dfop_sfo_saem) - vdiffr::expect_doppelganger("mixed model fit for saem object with saemix transformations", plot_biphasic_saem_s) - } - skip_on_travis() plot_biphasic_nlme <- function() plot(dfop_nlme_1) @@ -48,12 +43,6 @@ test_that("Plotting mkinfit, mmkin and mixed model objects is reproducible", { #plot_biphasic_mmkin <- function() plot(mixed(mmkin_biphasic)) # Biphasic fits with lots of data and fits have lots of potential for differences plot_biphasic_nlme <- function() plot(nlme_biphasic) - if (saemix_available) { - #plot_biphasic_saem_s <- function() plot(saem_biphasic_s) - plot_biphasic_saem_m <- function() plot(saem_biphasic_m) - - vdiffr::expect_doppelganger("mixed model fit for saem object with mkin transformations", plot_biphasic_saem_m) - } # different results when working with eigenvalues plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE) diff --git a/tests_slow.log b/tests_slow.log index 6337d80a..b5b646b6 100644 --- a/tests_slow.log +++ b/tests_slow.log @@ -1,11 +1,9 @@ +Loading required package: parallel ✔ | OK F W S | Context -✔ | 21 | Fitting of parent only models [39.2 s] -✔ | 15 | test_roundtrip_error_parameters [354.9 s] +✔ | 21 | Fitting of parent only models [24.2 s] +✔ | 15 | Roundtripping error model parameters [51.5 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 394.1 s +Duration: 75.7 s -OK: 36 -Failed: 0 -Warnings: 0 -Skipped: 0 +[ FAIL 0 | WARN 0 | SKIP 0 | PASS 36 ] diff --git a/vignettes/web_only/mkin_benchmarks.rda b/vignettes/web_only/mkin_benchmarks.rda index 4cd33ca1..b5f76c8e 100644 Binary files a/vignettes/web_only/mkin_benchmarks.rda and b/vignettes/web_only/mkin_benchmarks.rda differ -- cgit v1.2.1