From 804a5cb47fcdc823d41c585729ace151b283ca65 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 4 Aug 2023 13:27:41 +0200 Subject: Mesotrione data and vignette The vignette illustrates pH dependent degradation (covariate modelling) with some detail for the parent compound. Also, a bug in the saem method of the illparms function was fixed, which prevented to find ill-defined parameters in cases where e.g. slopes of covariate models have a negative estimate. --- DESCRIPTION | 4 +- NEWS.md | 8 + R/illparms.R | 4 +- inst/testdata/mesotrione_soil_efsa_2016.xlsx | Bin 0 -> 40810 bytes log/build.log | 2 +- log/check.log | 21 +- log/test.log | 54 +- .../_snaps/multistart/llhist-for-dfop-sfo-fit.svg | 43 +- ...t-for-saem-object-with-mkin-transformations.svg | 1086 ++++++++++---------- .../_snaps/multistart/parplot-for-dfop-sfo-fit.svg | 166 +-- .../plot/mixed-model-fit-for-nlme-object.svg | 20 +- tests/testthat/print_dfop_saem_1.txt | 10 +- vignettes/mesotrione_parent_2023_prebuilt.rnw | 7 + vignettes/prebuilt/2023_mesotrione_parent.pdf | Bin 0 -> 427383 bytes vignettes/prebuilt/2023_mesotrione_parent.rmd | 596 +++++++++++ 15 files changed, 1318 insertions(+), 703 deletions(-) create mode 100644 inst/testdata/mesotrione_soil_efsa_2016.xlsx create mode 100644 vignettes/mesotrione_parent_2023_prebuilt.rnw create mode 100644 vignettes/prebuilt/2023_mesotrione_parent.pdf create mode 100644 vignettes/prebuilt/2023_mesotrione_parent.rmd diff --git a/DESCRIPTION b/DESCRIPTION index a445a978..aaf1a04f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 1.2.4 -Date: 2023-05-16 +Version: 1.2.5 +Date: 2023-08-04 Authors@R: c( person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "johannes.ranke@jrwb.de", diff --git a/NEWS.md b/NEWS.md index b354afe6..fbf3da60 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# mkin 1.2.5 + +- 'vignettes/mesotrione_parent_2023.rnw': Prebuilt vignette showing how covariate modelling can be done for all relevant parent degradation models. + +- 'inst/testdata/mesotrione_soil_efsa_2016}.xlsx': Another example spreadsheets for use with 'read_spreadsheet()', featuring pH dependent degradation + +- R/illparms.R: Fix the detection of ill-defined slope or error model parameters for the case that the estimate is negative + # mkin 1.2.4 - R/endpoints.R: Fix the calculation of endpoints for user specified covariate values diff --git a/R/illparms.R b/R/illparms.R index eef4bd33..68a6bff6 100644 --- a/R/illparms.R +++ b/R/illparms.R @@ -106,7 +106,7 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod ill_parms <- c(ill_parms, names(which(ill_parms_random))) } if (errmod) { - ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0 + ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "upper"] > 0 ill_parms <- c(ill_parms, names(which(ill_parms_errmod))) } if (slopes) { @@ -115,7 +115,7 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod ci <- object$so@results@conf.int rownames(ci) <- ci$name slope_ci <- ci[slope_names, ] - ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "estimate"] > 0 + ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "upper"] > 0 ill_parms <- c(ill_parms, slope_names[ill_parms_slopes]) } } diff --git a/inst/testdata/mesotrione_soil_efsa_2016.xlsx b/inst/testdata/mesotrione_soil_efsa_2016.xlsx new file mode 100644 index 00000000..fa9ee7ae Binary files /dev/null and b/inst/testdata/mesotrione_soil_efsa_2016.xlsx differ diff --git a/log/build.log b/log/build.log index 30c5e432..717128fd 100644 --- a/log/build.log +++ b/log/build.log @@ -6,5 +6,5 @@ * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories Removed empty directory ‘mkin/vignettes/web_only’ -* building ‘mkin_1.2.4.tar.gz’ +* building ‘mkin_1.2.5.tar.gz’ diff --git a/log/check.log b/log/check.log index 970fbc7a..fedb0bd9 100644 --- a/log/check.log +++ b/log/check.log @@ -1,5 +1,5 @@ * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ -* using R version 4.3.0 Patched (2023-05-18 r84448) +* using R version 4.3.1 (2023-06-16) * using platform: x86_64-pc-linux-gnu (64-bit) * R was compiled by gcc (Debian 12.2.0-14) 12.2.0 @@ -9,9 +9,9 @@ * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘1.2.4’ +* this is package ‘mkin’ version ‘1.2.5’ * package encoding: UTF-8 -* checking CRAN incoming feasibility ... [2s/10s] Note_to_CRAN_maintainers +* checking CRAN incoming feasibility ... [3s/12s] Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke ’ * checking package namespace information ... OK * checking package dependencies ... OK @@ -45,7 +45,7 @@ Maintainer: ‘Johannes Ranke ’ * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK -* checking R code for possible problems ... OK +* checking R code for possible problems ... [15s/15s] OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -62,17 +62,22 @@ Maintainer: ‘Johannes Ranke ’ * checking sizes of PDF files under ‘inst/doc’ ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... OK +* checking examples ... [18s/18s] OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK * checking package vignettes in ‘inst/doc’ ... OK -* checking re-building of vignette outputs ... OK +* checking re-building of vignette outputs ... [17s/17s] OK * checking PDF version of manual ... OK -* checking HTML version of manual ... OK +* checking HTML version of manual ... NOTE +Skipping checking math rendering: package 'V8' unavailable * checking for non-standard things in the check directory ... OK * checking for detritus in the temp directory ... OK * DONE -Status: OK +Status: 1 NOTE +See + ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ +for details. + diff --git a/log/test.log b/log/test.log index b1b4ae31..2f1b6949 100644 --- a/log/test.log +++ b/log/test.log @@ -1,24 +1,24 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [1.5s] +✔ | 5 | Analytical solutions for coupled models [2.9s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE -✔ | 6 | Use of precompiled symbols in mkinpredict [3.1s] -✔ | 12 | Confidence intervals and p-values [0.4s] -✔ | 1 12 | Dimethenamid data from 2018 [12.5s] +✔ | 6 | Use of precompiled symbols in mkinpredict [5.7s] +✔ | 12 | Confidence intervals and p-values [1.0s] +✔ | 1 12 | Dimethenamid data from 2018 [33.1s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_dmta.R:88:3'): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ──────────────────────────────────────────────────────────────────────────────── -✔ | 14 | Error model fitting [2.3s] +✔ | 14 | Error model fitting [7.5s] ✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [0.3s] -✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.4s] -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.2s] -✔ | 1 | Fitting the logistic model [0.1s] -✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [18.9s] -✔ | 2 16 | Nonlinear mixed-effects models [148.4s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.5s] +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.7s] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] +✔ | 1 | Fitting the logistic model [0.2s] +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [43.6s] +✔ | 2 16 | Nonlinear mixed-effects models [364.4s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_mixed.R:80:3'): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve @@ -27,31 +27,31 @@ Skip ('test_mixed.R:133:3'): SFO-SFO saemix specific analytical solution work Reason: This is seldom used, so save some time ──────────────────────────────────────────────────────────────────────────────── ✔ | 3 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.3s] -✔ | 3 | mkinfit features [0.5s] -✔ | 8 | mkinmod model generation and printing -✔ | 3 | Model predictions with mkinpredict [0.1s] -✔ | 12 | Multistart method for saem.mmkin models [23.3s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.5s] -✔ | 9 | Nonlinear mixed-effects models with nlme [3.5s] -✔ | 15 | Plotting [4.5s] +✔ | 10 | Special cases of mkinfit calls [0.4s] +✔ | 3 | mkinfit features [0.7s] +✔ | 8 | mkinmod model generation and printing [0.2s] +✔ | 3 | Model predictions with mkinpredict [0.3s] +✔ | 12 | Multistart method for saem.mmkin models [80.2s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.0s] +✔ | 9 | Nonlinear mixed-effects models with nlme [9.1s] +✔ | 15 | Plotting [10.9s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [30.7s] +✔ | 1 36 | saemix parent models [77.0s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_saemix_parent.R:143:3'): We can also use mkin solution methods for saem Reason: This still takes almost 2.5 minutes although we do not solve ODEs ──────────────────────────────────────────────────────────────────────────────── -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [0.5s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.2s] ✔ | 11 | Processing of residue series -✔ | 10 | Fitting the SFORB model [1.7s] +✔ | 10 | Fitting the SFORB model [3.5s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 5 | Summary -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [0.8s] -✔ | 9 | Hypothesis tests [2.9s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [0.7s] +✔ | 5 | Summary [0.2s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.8s] +✔ | 9 | Hypothesis tests [6.7s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.0s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 259.6 s +Duration: 656.7 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) diff --git a/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg index 6015aed8..3b9d51fb 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg @@ -21,28 +21,26 @@ Frequency of log likelihoods - - - - + + + + + --1149.5 --1149.4 --1149.3 --1149.2 --1149.1 --1149.0 +-1149.30 +-1149.25 +-1149.20 +-1149.15 +-1149.10 +-1149.05 +-1149.00 - - 0 -1 -2 -3 -4 +1 +2 @@ -50,12 +48,13 @@ - - - - - - + + + + + + + original fit diff --git a/tests/testthat/_snaps/multistart/mixed-model-fit-for-saem-object-with-mkin-transformations.svg b/tests/testthat/_snaps/multistart/mixed-model-fit-for-saem-object-with-mkin-transformations.svg index 69fa6a4d..5e534dd1 100644 --- a/tests/testthat/_snaps/multistart/mixed-model-fit-for-saem-object-with-mkin-transformations.svg +++ b/tests/testthat/_snaps/multistart/mixed-model-fit-for-saem-object-with-mkin-transformations.svg @@ -249,7 +249,7 @@ - + @@ -342,7 +342,7 @@ - + @@ -856,38 +856,38 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -933,51 +933,51 @@ - - - + + + - - + + - - - - - - - - + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + @@ -1502,7 +1502,7 @@ - + @@ -1585,7 +1585,7 @@ - + @@ -1940,17 +1940,17 @@ - + - - - - + + + + 0 -10 -20 -30 -40 +10 +20 +30 +40 @@ -1983,290 +1983,290 @@ - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2274,47 +2274,47 @@ - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2323,174 +2323,174 @@ - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg index b01dac74..ed9168fb 100644 --- a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg +++ b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg @@ -25,104 +25,104 @@ - - - + + + - + - - - - - + + + + + - - - - - - + + + + + + - - - + + + - - - - - - + + + + + + - - - - - - - - - - + + + + + + + + + + - - - - + + + + - + - - - - - + + + + + - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - + - - - - - - + + + + + + - - + + - - - - - - - - - + + + + + + + + + - - + + diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg index 76fed0dc..c94012ce 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg @@ -813,7 +813,7 @@ - + @@ -833,7 +833,7 @@ - + @@ -900,7 +900,7 @@ - + @@ -921,7 +921,7 @@ - + @@ -964,8 +964,8 @@ - - + + @@ -1243,8 +1243,8 @@ - - + + @@ -1291,8 +1291,8 @@ - - + + diff --git a/tests/testthat/print_dfop_saem_1.txt b/tests/testthat/print_dfop_saem_1.txt index bdc40065..f7354320 100644 --- a/tests/testthat/print_dfop_saem_1.txt +++ b/tests/testthat/print_dfop_saem_1.txt @@ -13,11 +13,11 @@ Likelihood computed by importance sampling Fitted parameters: estimate lower upper -parent_0 99.92 98.77 101.06 -log_k1 -2.72 -2.95 -2.50 -log_k2 -4.14 -4.27 -4.01 -g_qlogis -0.35 -0.53 -0.16 -a.1 0.92 0.68 1.16 +parent_0 99.96 98.82 101.11 +log_k1 -2.71 -2.94 -2.49 +log_k2 -4.14 -4.26 -4.01 +g_qlogis -0.36 -0.54 -0.17 +a.1 0.93 0.69 1.17 b.1 0.05 0.04 0.06 SD.log_k1 0.37 0.23 0.51 SD.log_k2 0.23 0.14 0.31 diff --git a/vignettes/mesotrione_parent_2023_prebuilt.rnw b/vignettes/mesotrione_parent_2023_prebuilt.rnw new file mode 100644 index 00000000..f8a7a097 --- /dev/null +++ b/vignettes/mesotrione_parent_2023_prebuilt.rnw @@ -0,0 +1,7 @@ +\documentclass{article} +\usepackage{pdfpages} +%\VignetteIndexEntry{Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione} + +\begin{document} +\includepdf[pages=-, fitpaper=true]{prebuilt/2023_mesotrione_parent.pdf} +\end{document} diff --git a/vignettes/prebuilt/2023_mesotrione_parent.pdf b/vignettes/prebuilt/2023_mesotrione_parent.pdf new file mode 100644 index 00000000..47d18761 Binary files /dev/null and b/vignettes/prebuilt/2023_mesotrione_parent.pdf differ diff --git a/vignettes/prebuilt/2023_mesotrione_parent.rmd b/vignettes/prebuilt/2023_mesotrione_parent.rmd new file mode 100644 index 00000000..6ab126d9 --- /dev/null +++ b/vignettes/prebuilt/2023_mesotrione_parent.rmd @@ -0,0 +1,596 @@ +--- +title: "Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione" +author: Johannes Ranke +date: Last change on 4 August 2023, last compiled on `r format(Sys.time(), "%e + %B %Y")` +output: + pdf_document: + extra_dependencies: ["float", "listing"] +toc: yes +geometry: margin=2cm +--- + +```{r setup, echo = FALSE, cache = FALSE} +options(width = 80) # For summary listings +knitr::opts_chunk$set( + cache = TRUE, + comment = "", tidy = FALSE, + fig.pos = "H", fig.align = "center" +) +``` + +\clearpage + +# Introduction + +The purpose of this document is to test demonstrate how nonlinear hierarchical +models (NLHM) based on the parent degradation models SFO, FOMC, DFOP and HS +can be fitted with the mkin package, also considering the influence of +covariates like soil pH on different degradation parameters. Because in some +other case studies, the SFORB parameterisation of biexponential decline has +shown some advantages over the DFOP parameterisation, SFORB was included +in the list of tested models as well. + +The mkin package is used in version `r packageVersion("mkin")`, which is +contains the functions that were used for the evaluations. The `saemix` package +is used as a backend for fitting the NLHM, but is also loaded to make the +convergence plot function available. + +This document is processed with the `knitr` package, which also provides the +`kable` function that is used to improve the display of tabular data in R +markdown documents. For parallel processing, the `parallel` package is used. + +```{r packages, cache = FALSE, message = FALSE} +library(mkin) +library(knitr) +library(saemix) +library(parallel) +n_cores <- detectCores() +if (Sys.info()["sysname"] == "Windows") { + cl <- makePSOCKcluster(n_cores) +} else { + cl <- makeForkCluster(n_cores) +} +``` + +\clearpage + +## Test data + +```{r data} +data_file <- system.file( + "testdata", "mesotrione_soil_efsa_2016.xlsx", package = "mkin") +meso_ds <- read_spreadsheet(data_file, parent_only = TRUE) +``` + +The following tables show the covariate data and the `r length(meso_ds)` +datasets that were read in from the spreadsheet file. + +```{r show-covar-data, dependson = "data", results = "asis"} +pH <- attr(meso_ds, "covariates") +kable(pH, caption = "Covariate data") +``` + +\clearpage + +```{r show-data, dependson = "data", results = "asis"} +for (ds_name in names(meso_ds)) { + print( + kable(mkin_long_to_wide(meso_ds[[ds_name]]), + caption = paste("Dataset", ds_name), + booktabs = TRUE, row.names = FALSE)) +} +``` + +\clearpage + +# Separate evaluations + +In order to obtain suitable starting parameters for the NLHM fits, +separate fits of the five models to the data for each +soil are generated using the `mmkin` function from the mkin package. +In a first step, constant variance is assumed. Convergence +is checked with the `status` function. + +```{r f-sep-const, dependson = "data"} +deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB", "HS") +f_sep_const <- mmkin( + deg_mods, + meso_ds, + error_model = "const", + cluster = cl, + quiet = TRUE) +``` + +```{r dependson = "f-sep-const"} +status(f_sep_const[, 1:5]) |> kable() +status(f_sep_const[, 6:18]) |> kable() +``` + +In the tables above, OK indicates convergence and C indicates failure to +converge. Most separate fits with constant variance converged, with the +exception of two FOMC fits, one SFORB fit and one HS fit. + +```{r f-sep-tc, dependson = "f-sep-const"} +f_sep_tc <- update(f_sep_const, error_model = "tc") +``` + +```{r dependson = "f-sep-tc"} +status(f_sep_tc[, 1:5]) |> kable() +status(f_sep_tc[, 6:18]) |> kable() +``` + +With the two-component error model, the set of fits that did not converge +is larger, with convergence problems appearing for a number of non-SFO fits. + +\clearpage + +# Hierarchical model fits without covariate effect + +The following code fits hierarchical kinetic models for the ten combinations of +the five different degradation models with the two different error models in +parallel. + +```{r f-saem-1, dependson = c("f-sep-const", "f-sep-tc")} +f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cluster = cl) +status(f_saem_1) |> kable() +``` + +All fits terminate without errors (status OK). + +```{r dependson = "f-saem-1"} +anova(f_saem_1) |> kable(digits = 1) +``` +The model comparisons show that the fits with constant variance are consistently +preferable to the corresponding fits with two-component error for these data. +This is confirmed by the fact that the parameter `b.1` (the relative standard +deviation in the fits obtained with the saemix package), is ill-defined in all +fits. + +```{r dependson = "f-saem-1"} +illparms(f_saem_1) |> kable() +``` + +For obtaining fits with only well-defined random effects, we update +the set of fits, excluding random effects that were ill-defined +according to the `illparms` function. + +```{r f-saem-2, dependson = "f-saem-1"} +f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1)) +status(f_saem_2) |> kable() +``` + +The updated fits terminate without errors. + +```{r dependson = "f-saem-2"} +illparms(f_saem_2) |> kable() +``` + +No ill-defined errors remain in the fits with constant variance. + +\clearpage + +# Hierarchical model fits with covariate effect + +In the following sections, hierarchical fits including a model for the +influence of pH on selected degradation parameters are shown for all parent +models. Constant variance is selected as the error model based on the fits +without covariate effects. Random effects that were ill-defined in the fits +without pH influence are excluded. A potential influence of the soil pH is only +included for parameters with a well-defined random effect, because experience +has shown that for such parameters no significant pH effect could be found. + +## SFO + +```{r sfo-pH, dependson = "f-sep-const"} +sfo_pH <- saem(f_sep_const["SFO", ], no_random_effect = "meso_0", covariates = pH, + covariate_models = list(log_k_meso ~ pH)) +``` + +```{r dependson = "sfo-pH"} +summary(sfo_pH)$confint_trans |> kable(digits = 2) +``` + +The parameter showing the pH influence in the above table is +`beta_pH(log_k_meso)`. Its confidence interval does not include zero, +indicating that the influence of soil pH on the log of the degradation rate +constant is significantly greater than zero. + + +```{r dependson = "sfo-pH"} +anova(f_saem_2[["SFO", "const"]], sfo_pH, test = TRUE) +``` + +The comparison with the SFO fit without covariate effect confirms that +considering the soil pH improves the model, both by comparison of AIC and BIC +and by the likelihood ratio test. + +\clearpage + +```{r dependson = "sfo-pH", plot.height = 5} +plot(sfo_pH) +``` + +Endpoints for a model with covariates are by default calculated for +the median of the covariate values. This quantile can be adapted, +or a specific covariate value can be given as shown below. + +```{r dependson = "sfo-pH", plot.height = 5} +endpoints(sfo_pH) +endpoints(sfo_pH, covariate_quantile = 0.9) +endpoints(sfo_pH, covariates = c(pH = 7.0)) +``` + +\clearpage + +## FOMC + +```{r fomc-pH, dependson = "f-sep-const"} +fomc_pH <- saem(f_sep_const["FOMC", ], no_random_effect = "meso_0", covariates = pH, + covariate_models = list(log_alpha ~ pH)) +``` + +```{r dependson = "fomc-pH"} +summary(fomc_pH)$confint_trans |> kable(digits = 2) +``` + +As in the case of SFO, the confidence interval of the slope parameter +(here `beta_pH(log_alpha)`) quantifying the influence of soil pH +does not include zero, and the model comparison clearly indicates +that the model with covariate influence is preferable. +However, the random effect for `alpha` is not well-defined any +more after inclusion of the covariate effect (the confidence +interval of `SD.log_alpha` includes zero). + +```{r dependson = "fomc-pH"} +illparms(fomc_pH) +``` + +Therefore, the model is updated without this random effect, and +no ill-defined parameters remain. + +```{r fomc-pH-2, dependson = "fomc_pH"} +fomc_pH_2 <- update(fomc_pH, no_random_effect = c("meso_0", "log_alpha")) +illparms(fomc_pH_2) +``` + +```{r dependson = "fomc-pH-2"} +anova(f_saem_2[["FOMC", "const"]], fomc_pH, fomc_pH_2, test = TRUE) +``` + +Model comparison indicates that including pH dependence significantly improves +the fit, and that the reduced model with covariate influence results in +the most preferable FOMC fit. + +```{r dependson = "fomc-pH"} +summary(fomc_pH_2)$confint_trans |> kable(digits = 2) +``` + +\clearpage + +```{r dependson = "fomc-pH", plot.height = 5} +plot(fomc_pH_2) +``` + +```{r dependson = "fomc-pH", plot.height = 5} +endpoints(fomc_pH_2) +endpoints(fomc_pH_2, covariates = c(pH = 7)) +``` + +\clearpage + +## DFOP + +In the DFOP fits without covariate effects, random effects for two degradation +parameters (`k2` and `g`) were identifiable. + +```{r dependson = "f-saem-2"} +summary(f_saem_2[["DFOP", "const"]])$confint_trans |> kable(digits = 2) +``` + +A fit with pH dependent degradation parameters was obtained by excluding +the same random effects as in the refined DFOP fit without covariate influence, +and including covariate models for the two identifiable parameters `k2` and `g`. + +```{r dfop-pH, dependson = "f-sep-const"} +dfop_pH <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"), + covariates = pH, + covariate_models = list(log_k2 ~ pH, g_qlogis ~ pH)) +``` + +The corresponding parameters for +the influence of soil pH are `beta_pH(log_k2)` for the influence of soil pH on +`k2`, and `beta_pH(g_qlogis)` for its influence on `g`. + +```{r dependson = "dfop-pH"} +summary(dfop_pH)$confint_trans |> kable(digits = 2) +illparms(dfop_pH) +``` + +Confidence intervals for neither of them include zero, indicating a significant +difference from zero. However, the random effect for `g` is now ill-defined. +The fit is updated without this ill-defined random effect. + +```{r dfop-pH-2, dependson = "dfop-pH"} +dfop_pH_2 <- update(dfop_pH, + no_random_effect = c("meso_0", "log_k1", "g_qlogis")) +illparms(dfop_pH_2) +``` + +Now, the slope parameter for the pH effect on `g` is ill-defined. +Therefore, another attempt is made without the corresponding covariate model. + +```{r dfop-pH-3, dependson = "f-sep-const"} +dfop_pH_3 <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"), + covariates = pH, + covariate_models = list(log_k2 ~ pH)) +illparms(dfop_pH_3) +``` +As the random effect for `g` is again ill-defined, the fit is repeated without it. + +```{r dfop-pH-4, dependson = "dfop-pH-3"} +dfop_pH_4 <- update(dfop_pH_3, no_random_effect = c("meso_0", "log_k1", "g_qlogis")) +illparms(dfop_pH_4) +``` + +While no ill-defined parameters remain, model comparison suggests that the previous +model `dfop_pH_2` with two pH dependent parameters is preferable, based on +information criteria as well as based on the likelihood ratio test. + +```{r dependson = "dfop-pH-4"} +anova(f_saem_2[["DFOP", "const"]], dfop_pH, dfop_pH_2, dfop_pH_3, dfop_pH_4) +anova(dfop_pH_2, dfop_pH_4, test = TRUE) +``` + +When focussing on parameter identifiability using the test if the confidence +interval includes zero, `dfop_pH_4` would still be the preferred model. +However, it should be kept in mind that parameter confidence intervals are +constructed using a simple linearisation of the likelihood. As the confidence +interval of the random effect for `g` only marginally includes zero, +it is suggested that this is acceptable, and that `dfop_pH_2` can be considered +the most preferable model. + +\clearpage + +```{r dependson = "dfop-pH-2", plot.height = 5} +plot(dfop_pH_2) +``` + +```{r dependson = "dfop-pH-2", plot.height = 5} +endpoints(dfop_pH_2) +endpoints(dfop_pH_2, covariates = c(pH = 7)) +``` + +\clearpage + +## SFORB + +```{r sforb-pH, dependson = "f-sep-const"} +sforb_pH <- saem(f_sep_const["SFORB", ], no_random_effect = c("meso_free_0", "log_k_meso_free_bound"), + covariates = pH, + covariate_models = list(log_k_meso_free ~ pH, log_k_meso_bound_free ~ pH)) +``` + +```{r dependson = "sforb-pH"} +summary(sforb_pH)$confint_trans |> kable(digits = 2) +``` + +The confidence interval of `beta_pH(log_k_meso_bound_free)` includes zero, +indicating that the influence of soil pH on `k_meso_bound_free` cannot reliably +be quantified. Also, the confidence interval for the random effect on this +parameter (`SD.log_k_meso_bound_free`) includes zero. + +Using the `illparms` function, these ill-defined parameters can be found +more conveniently. + +```{r dependson = "sforb-pH"} +illparms(sforb_pH) +``` + +To remove the ill-defined parameters, a second variant of the SFORB model +with pH influence is fitted. No ill-defined parameters remain. + +```{r sforb-pH-2, dependson = "f-sforb-pH"} +sforb_pH_2 <- update(sforb_pH, + no_random_effect = c("meso_free_0", "log_k_meso_free_bound", "log_k_meso_bound_free"), + covariate_models = list(log_k_meso_free ~ pH)) +illparms(sforb_pH_2) +``` + +The model comparison of the SFORB fits includes the refined model without +covariate effect, and both versions of the SFORB fit with covariate effect. + +```{r dependson = "sforb-pH-2"} +anova(f_saem_2[["SFORB", "const"]], sforb_pH, sforb_pH_2, test = TRUE) +``` +The first model including pH influence is preferable based on information criteria +and the likelihood ratio test. However, as it is not fully identifiable, +the second model is selected. + +```{r dependson = "sforb-pH-2"} +summary(sforb_pH_2)$confint_trans |> kable(digits = 2) +``` +\clearpage + +```{r dependson = "sforb-pH-2", plot.height = 5} +plot(sforb_pH_2) +``` + +```{r dependson = "sforb-pH-2", plot.height = 5} +endpoints(sforb_pH_2) +endpoints(sforb_pH_2, covariates = c(pH = 7)) +``` + +\clearpage + +## HS + +```{r hs-pH, dependson = "f-sep-const"} +hs_pH <- saem(f_sep_const["HS", ], no_random_effect = c("meso_0"), + covariates = pH, + covariate_models = list(log_k1 ~ pH, log_k2 ~ pH, log_tb ~ pH)) +``` + +```{r dependson = "hs-pH"} +summary(hs_pH)$confint_trans |> kable(digits = 2) +illparms(hs_pH) +``` + +According to the output of the `illparms` function, the random effect on +the break time `tb` cannot reliably be quantified, neither can the influence of +soil pH on `tb`. The fit is repeated without the corresponding covariate +model, and no ill-defined parameters remain. + +```{r hs-pH-2, dependson = "hs-pH"} +hs_pH_2 <- update(hs_pH, covariate_models = list(log_k1 ~ pH, log_k2 ~ pH)) +illparms(hs_pH_2) +``` + +Model comparison confirms that this model is preferable to the fit without +covariate influence, and also to the first version with covariate influence. + +```{r dependson = c("hs-pH-2", "hs-pH-3")} +anova(f_saem_2[["HS", "const"]], hs_pH, hs_pH_2, test = TRUE) +``` + +```{r dependson = "hs-pH-2"} +summary(hs_pH_2)$confint_trans |> kable(digits = 2) +``` + +\clearpage + +```{r dependson = "hs-pH-2", plot.height = 5} +plot(hs_pH_2) +``` + +```{r dependson = "hs-pH-2", plot.height = 5} +endpoints(hs_pH_2) +endpoints(hs_pH_2, covariates = c(pH = 7)) +``` + +\clearpage + +## Comparison across parent models + +After model reduction for all models with pH influence, they are compared with +each other. + +```{r, dependson = c("sfo-pH-2", "fomc-pH-2", "dfop-pH-4", "sforb-pH-1", "hs-pH-3")} +anova(sfo_pH, fomc_pH_2, dfop_pH_2, dfop_pH_4, sforb_pH_2, hs_pH_2) +``` + +The DFOP model with pH influence on `k2` and `g` and a random effect only on +`k2` is finally selected as the best fit. + +The endpoints resulting from this model are listed below. Please refer to the +Appendix for a detailed listing. + +```{r, dependson = "dfop-pH-2"} +endpoints(dfop_pH_2) +endpoints(dfop_pH_2, covariates = c(pH = 7)) +``` + +# Conclusions + +These evaluations demonstrate that covariate effects can be included +for all types of parent degradation models. These models can then +be further refined to make them fully identifiable. + +\clearpage + +# Appendix + +## Hierarchical fit listings + +### Fits without covariate effects + +```{r, cache = FALSE, results = "asis", echo = FALSE} +errmods <- c(const = "constant variance", tc = "two-component error") +for (deg_mod in deg_mods) { + for (err_mod in c("const")) { + fit <- f_saem_1[[deg_mod, err_mod]] + if (!inherits(fit$so, "try-error")) { + caption <- paste("Hierarchical", deg_mod, "fit with", errmods[err_mod]) + tex_listing(fit, caption) + } + } +} +``` + +### Fits with covariate effects + +```{r listing-sfo, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(sfo_pH, "Hierarchichal SFO fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(fomc_pH, "Hierarchichal FOMC fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(fomc_pH_2, "Refined hierarchichal FOMC fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(dfop_pH, "Hierarchichal DFOP fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(dfop_pH_2, "Refined hierarchical DFOP fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(dfop_pH_4, "Further refined hierarchical DFOP fit with pH influence") +``` + +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(sforb_pH, "Hierarchichal SFORB fit with pH influence") +``` +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(sforb_pH_2, "Refined hierarchichal SFORB fit with pH influence") +``` +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(hs_pH, "Hierarchichal HS fit with pH influence") +``` +\clearpage + +```{r, cache = FALSE, results = "asis", echo = FALSE} +tex_listing(hs_pH_2, "Refined hierarchichal HS fit with pH influence") +``` + +\clearpage + +## Session info + +```{r, echo = FALSE, cache = FALSE} +parallel::stopCluster(cl = cl) +sessionInfo() +``` + + +## Hardware info + +```{r, echo = FALSE} +if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) { + cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]])) +} +if(!inherits(try(meminfo <- readLines("/proc/meminfo")), "try-error")) { + cat(gsub("model name\t: ", "System memory: ", meminfo[grep("MemTotal", meminfo)[1]])) +} +``` -- cgit v1.2.1 From b1754e3d1112ea04e2ee4c393581362fb8f605f7 Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Mon, 7 Aug 2023 11:06:17 +0200 Subject: Typos --- R/read_spreadsheet.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/read_spreadsheet.R b/R/read_spreadsheet.R index 1c5508da..bad93d43 100644 --- a/R/read_spreadsheet.R +++ b/R/read_spreadsheet.R @@ -12,7 +12,7 @@ #' to contain name and acronym of the parent compound. #' #' The dataset sheets should be named using the dataset numbers read in from -#' the 'Datasets' sheet, i.e. '1', '2', ... . In each dataset sheet, name +#' the 'Datasets' sheet, i.e. '1', '2', ... . In each dataset sheet, the name #' of the observed variable (e.g. the acronym of the parent compound or #' one of its transformation products) should be in the first column, #' the time values should be in the second colum, and the observed value @@ -26,7 +26,7 @@ #' Their names should preferably not contain special characters like spaces, #' so they can be easily used for specifying covariate models. #' -#' An similar data structure is defined as the R6 class [mkindsg], but +#' A similar data structure is defined as the R6 class [mkindsg], but #' is probably more complicated to use. #' #' @param path Absolute or relative path to the spreadsheet file -- cgit v1.2.1 From 8155e46f24ad2aa7562ced1c52333fa95d072d14 Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Tue, 8 Aug 2023 15:20:13 +0200 Subject: RStudio project infrastructure --- .Rbuildignore | 2 ++ .gitignore | 1 + mkin.Rproj | 22 ++++++++++++++++++++++ 3 files changed, 25 insertions(+) create mode 100644 mkin.Rproj diff --git a/.Rbuildignore b/.Rbuildignore index 0ee54702..e7cbe786 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -33,3 +33,5 @@ ^vignettes/FOCUS_Z.tex$ ^vignettes/mkin.tex$ ^vignettes/web_only/.*$ +^.*\.Rproj$ +^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore index babe89e7..b383dc63 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,4 @@ vignettes/*_files/ vignettes/prebuilt/*_cache vignettes/prebuilt/*_files vignettes/prebuilt/*_dlls +.Rproj.user diff --git a/mkin.Rproj b/mkin.Rproj new file mode 100644 index 00000000..fac143c3 --- /dev/null +++ b/mkin.Rproj @@ -0,0 +1,22 @@ +Version: 1.0 + +RestoreWorkspace: No +SaveWorkspace: No +AlwaysSaveHistory: No + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: knitr +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Makefile + +UseNativePipeOperator: Yes + +SpellingDictionary: en_GB -- cgit v1.2.1 From ec04610fca09cedad7bc8c367ae6884a1bd7b7f6 Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Tue, 8 Aug 2023 15:39:08 +0200 Subject: Fix a logical error in the mesotrione vignette --- vignettes/prebuilt/2023_mesotrione_parent.pdf | Bin 427383 -> 434484 bytes vignettes/prebuilt/2023_mesotrione_parent.rmd | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/prebuilt/2023_mesotrione_parent.pdf b/vignettes/prebuilt/2023_mesotrione_parent.pdf index 47d18761..a7ecfe3c 100644 Binary files a/vignettes/prebuilt/2023_mesotrione_parent.pdf and b/vignettes/prebuilt/2023_mesotrione_parent.pdf differ diff --git a/vignettes/prebuilt/2023_mesotrione_parent.rmd b/vignettes/prebuilt/2023_mesotrione_parent.rmd index 6ab126d9..f11effc7 100644 --- a/vignettes/prebuilt/2023_mesotrione_parent.rmd +++ b/vignettes/prebuilt/2023_mesotrione_parent.rmd @@ -178,7 +178,7 @@ models. Constant variance is selected as the error model based on the fits without covariate effects. Random effects that were ill-defined in the fits without pH influence are excluded. A potential influence of the soil pH is only included for parameters with a well-defined random effect, because experience -has shown that for such parameters no significant pH effect could be found. +has shown that only for such parameters a significant pH effect could be found. ## SFO -- cgit v1.2.1 From 71d2f67189eaaad9b8d61d0f5c709f8d7be65a68 Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Tue, 8 Aug 2023 16:56:59 +0200 Subject: Update docs, check and test on RStudio server --- log/build.log | 2 ++ log/check.log | 20 ++++++------- log/test.log | 76 +++++++++++++++++++++---------------------------- man/read_spreadsheet.Rd | 4 +-- mkin.Rproj | 3 ++ 5 files changed, 49 insertions(+), 56 deletions(-) diff --git a/log/build.log b/log/build.log index 717128fd..11086361 100644 --- a/log/build.log +++ b/log/build.log @@ -7,4 +7,6 @@ * checking for empty or unneeded directories Removed empty directory ‘mkin/vignettes/web_only’ * building ‘mkin_1.2.5.tar.gz’ +Warning: invalid uid value replaced by that for user 'nobody' +Warning: invalid gid value replaced by that for user 'nobody' diff --git a/log/check.log b/log/check.log index fedb0bd9..64457e57 100644 --- a/log/check.log +++ b/log/check.log @@ -1,17 +1,17 @@ -* using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ +* using log directory ‘/home/agsad.admin.ch/f80868656/projects/mkin/mkin.Rcheck’ * using R version 4.3.1 (2023-06-16) * using platform: x86_64-pc-linux-gnu (64-bit) * R was compiled by - gcc (Debian 12.2.0-14) 12.2.0 - GNU Fortran (Debian 12.2.0-14) 12.2.0 -* running under: Debian GNU/Linux 12 (bookworm) + gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0 + GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0 +* running under: Ubuntu 22.04.3 LTS * using session charset: UTF-8 * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package * this is package ‘mkin’ version ‘1.2.5’ * package encoding: UTF-8 -* checking CRAN incoming feasibility ... [3s/12s] Note_to_CRAN_maintainers +* checking CRAN incoming feasibility ... [5s/26s] Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke ’ * checking package namespace information ... OK * checking package dependencies ... OK @@ -45,7 +45,7 @@ Maintainer: ‘Johannes Ranke ’ * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK -* checking R code for possible problems ... [15s/15s] OK +* checking R code for possible problems ... [14s/15s] OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -62,22 +62,22 @@ Maintainer: ‘Johannes Ranke ’ * checking sizes of PDF files under ‘inst/doc’ ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... [18s/18s] OK +* checking examples ... [20s/22s] OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK * checking package vignettes in ‘inst/doc’ ... OK -* checking re-building of vignette outputs ... [17s/17s] OK +* checking re-building of vignette outputs ... [17s/18s] OK * checking PDF version of manual ... OK * checking HTML version of manual ... NOTE -Skipping checking math rendering: package 'V8' unavailable +Skipping checking HTML validation: no command 'tidy' found * checking for non-standard things in the check directory ... OK * checking for detritus in the temp directory ... OK * DONE Status: 1 NOTE See - ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ + ‘/home/agsad.admin.ch/f80868656/projects/mkin/mkin.Rcheck/00check.log’ for details. diff --git a/log/test.log b/log/test.log index 2f1b6949..3bca2e6a 100644 --- a/log/test.log +++ b/log/test.log @@ -4,59 +4,47 @@ ✔ | 5 | Analytical solutions for coupled models [2.9s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE -✔ | 6 | Use of precompiled symbols in mkinpredict [5.7s] -✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [33.1s] -──────────────────────────────────────────────────────────────────────────────── -Skip ('test_dmta.R:88:3'): Different backends get consistent results for SFO-SFO3+, dimethenamid data -Reason: Fitting this ODE model with saemix takes about 15 minutes on my system -──────────────────────────────────────────────────────────────────────────────── -✔ | 14 | Error model fitting [7.5s] +✔ | 6 | Use of precompiled symbols in mkinpredict [6.6s] +✔ | 12 | Confidence intervals and p-values [1.1s] +✔ | 1 12 | Dimethenamid data from 2018 [32.2s] +✔ | 14 | Error model fitting [5.4s] ✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [0.5s] -✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.7s] -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] -✔ | 1 | Fitting the logistic model [0.2s] -✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [43.6s] -✔ | 2 16 | Nonlinear mixed-effects models [364.4s] -──────────────────────────────────────────────────────────────────────────────── -Skip ('test_mixed.R:80:3'): saemix results are reproducible for biphasic fits -Reason: Fitting with saemix takes around 10 minutes when using deSolve - -Skip ('test_mixed.R:133:3'): SFO-SFO saemix specific analytical solution work -Reason: This is seldom used, so save some time -──────────────────────────────────────────────────────────────────────────────── +✔ | 4 | Calculation of FOCUS chi2 error levels +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) +✔ | 4 | Test fitting the decline of metabolites from their maximum +✔ | 1 | Fitting the logistic model +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [46.9s] +✔ | 2 16 | Nonlinear mixed-effects models [399.6s] ✔ | 3 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.4s] -✔ | 3 | mkinfit features [0.7s] -✔ | 8 | mkinmod model generation and printing [0.2s] -✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 12 | Multistart method for saem.mmkin models [80.2s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.0s] -✔ | 9 | Nonlinear mixed-effects models with nlme [9.1s] +✔ | 10 | Special cases of mkinfit calls +✔ | 3 | mkinfit features +✔ | 8 | mkinmod model generation and printing +✔ | 3 | Model predictions with mkinpredict +✔ | 12 | Multistart method for saem.mmkin models [50.1s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.8s] +✔ | 9 | Nonlinear mixed-effects models with nlme [9.3s] ✔ | 15 | Plotting [10.9s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [77.0s] -──────────────────────────────────────────────────────────────────────────────── -Skip ('test_saemix_parent.R:143:3'): We can also use mkin solution methods for saem -Reason: This still takes almost 2.5 minutes although we do not solve ODEs -──────────────────────────────────────────────────────────────────────────────── +✔ | 1 36 | saemix parent models [68.4s] ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.2s] ✔ | 11 | Processing of residue series -✔ | 10 | Fitting the SFORB model [3.5s] +✔ | 10 | Fitting the SFORB model [3.8s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 5 | Summary [0.2s] -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.8s] -✔ | 9 | Hypothesis tests [6.7s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.0s] +✔ | 5 | Summary +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.9s] +✔ | 9 | Hypothesis tests [6.4s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [1.4s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 656.7 s +Duration: 655.6 s -── Skipped tests ────────────────────────────────────────────────────────────── -• Fitting this ODE model with saemix takes about 15 minutes on my system (1) -• Fitting with saemix takes around 10 minutes when using deSolve (1) -• This is seldom used, so save some time (1) -• This still takes almost 2.5 minutes although we do not solve ODEs (1) +── Skipped tests (4) ─────────────────────────────────────────────────────────── +• Fitting this ODE model with saemix takes about 15 minutes on my system (1): + 'test_dmta.R:88:3' +• Fitting with saemix takes around 10 minutes when using deSolve (1): + 'test_mixed.R:80:3' +• This is seldom used, so save some time (1): 'test_mixed.R:133:3' +• This still takes almost 2.5 minutes although we do not solve ODEs (1): + 'test_saemix_parent.R:143:3' [ FAIL 0 | WARN 0 | SKIP 4 | PASS 281 ] diff --git a/man/read_spreadsheet.Rd b/man/read_spreadsheet.Rd index 41c32108..eefda692 100644 --- a/man/read_spreadsheet.Rd +++ b/man/read_spreadsheet.Rd @@ -36,7 +36,7 @@ The first row read after the header read in from this sheet is assumed to contain name and acronym of the parent compound. The dataset sheets should be named using the dataset numbers read in from -the 'Datasets' sheet, i.e. '1', '2', ... . In each dataset sheet, name +the 'Datasets' sheet, i.e. '1', '2', ... . In each dataset sheet, the name of the observed variable (e.g. the acronym of the parent compound or one of its transformation products) should be in the first column, the time values should be in the second colum, and the observed value @@ -50,6 +50,6 @@ column and the column must have the same name as the second column in Their names should preferably not contain special characters like spaces, so they can be easily used for specifying covariate models. -An similar data structure is defined as the R6 class \link{mkindsg}, but +A similar data structure is defined as the R6 class \link{mkindsg}, but is probably more complicated to use. } diff --git a/mkin.Rproj b/mkin.Rproj index fac143c3..beb8549b 100644 --- a/mkin.Rproj +++ b/mkin.Rproj @@ -19,4 +19,7 @@ BuildType: Makefile UseNativePipeOperator: Yes +MarkdownWrap: Column +MarkdownWrapAtColumn: 80 + SpellingDictionary: en_GB -- cgit v1.2.1 From bd13f94b35e325da350a247c88b249dba62fe22f Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Wed, 9 Aug 2023 16:03:44 +0200 Subject: Don't install saemix 3.2 from github any more Emmanuelle has made an official release to CRAN in June --- .gitignore | 1 + .travis.yml | 1 - DESCRIPTION | 2 +- log/test.log | 34 +++++++++++++++++----------------- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/.gitignore b/.gitignore index b383dc63..4a79508d 100644 --- a/.gitignore +++ b/.gitignore @@ -33,3 +33,4 @@ vignettes/prebuilt/*_cache vignettes/prebuilt/*_files vignettes/prebuilt/*_dlls .Rproj.user +.Rhistory diff --git a/.travis.yml b/.travis.yml index ec904005..f6beb25b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,7 +19,6 @@ before_install: install: - ./run.sh install_all - - ./run.sh install_github jranke/saemixextension@installable_dev_version script: - travis_wait 30 ./run.sh run_tests diff --git a/DESCRIPTION b/DESCRIPTION index aaf1a04f..ea5ee18d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data Version: 1.2.5 -Date: 2023-08-04 +Date: 2023-08-09 Authors@R: c( person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "johannes.ranke@jrwb.de", diff --git a/log/test.log b/log/test.log index 3bca2e6a..c872cb1d 100644 --- a/log/test.log +++ b/log/test.log @@ -1,42 +1,42 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [2.9s] +✔ | 5 | Analytical solutions for coupled models [3.1s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE -✔ | 6 | Use of precompiled symbols in mkinpredict [6.6s] -✔ | 12 | Confidence intervals and p-values [1.1s] -✔ | 1 12 | Dimethenamid data from 2018 [32.2s] -✔ | 14 | Error model fitting [5.4s] +✔ | 6 | Use of precompiled symbols in mkinpredict [5.3s] +✔ | 12 | Confidence intervals and p-values +✔ | 1 12 | Dimethenamid data from 2018 [31.5s] +✔ | 14 | Error model fitting [5.6s] ✔ | 5 | Time step normalisation ✔ | 4 | Calculation of FOCUS chi2 error levels ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ✔ | 4 | Test fitting the decline of metabolites from their maximum ✔ | 1 | Fitting the logistic model -✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [46.9s] -✔ | 2 16 | Nonlinear mixed-effects models [399.6s] +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [43.1s] +✔ | 2 16 | Nonlinear mixed-effects models [441.3s] ✔ | 3 | Test dataset classes mkinds and mkindsg ✔ | 10 | Special cases of mkinfit calls ✔ | 3 | mkinfit features ✔ | 8 | mkinmod model generation and printing ✔ | 3 | Model predictions with mkinpredict -✔ | 12 | Multistart method for saem.mmkin models [50.1s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.8s] -✔ | 9 | Nonlinear mixed-effects models with nlme [9.3s] -✔ | 15 | Plotting [10.9s] +✔ | 12 | Multistart method for saem.mmkin models [71.6s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.7s] +✔ | 9 | Nonlinear mixed-effects models with nlme [10.1s] +✔ | 15 | Plotting [12.1s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [68.4s] -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.2s] +✔ | 1 36 | saemix parent models [81.8s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s] ✔ | 11 | Processing of residue series ✔ | 10 | Fitting the SFORB model [3.8s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.9s] -✔ | 9 | Hypothesis tests [6.4s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [1.4s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s] +✔ | 9 | Hypothesis tests [7.3s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [1.5s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 655.6 s +Duration: 731.0 s ── Skipped tests (4) ─────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1): -- cgit v1.2.1 From 60a8843ac44b97ef819d30968564935375ca470c Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Wed, 9 Aug 2023 16:09:50 +0200 Subject: Depend on saemix 3.2 to control test results --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index ea5ee18d..916d1cdc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,7 +24,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, particular purpose. Depends: R (>= 2.15.1), Imports: stats, graphics, methods, parallel, deSolve (>= 1.35), R6, inline (>= 0.3.19), - numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), saemix (>= 3.1), rlang, vctrs + numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), saemix (>= 3.2), rlang, vctrs Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, benchmarkme, tibble, stats4, readxl License: GPL -- cgit v1.2.1 From 9f76ca630be7043333b9e1cc51e0062d87cf76d2 Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Wed, 9 Aug 2023 16:29:58 +0200 Subject: Try jammy on travis --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index f6beb25b..d32e9dde 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,6 @@ language: c sudo: required -dist: focal +dist: jammy jobs: include: -- cgit v1.2.1 From 1bda49a550ad5ed66a4eade350850e33f1033a7a Mon Sep 17 00:00:00 2001 From: Ranke Johannes Date: Wed, 9 Aug 2023 16:40:45 +0200 Subject: Go back to focal, skip drat repo on travis --- .travis.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index d32e9dde..683ca4ee 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,6 @@ language: c sudo: required -dist: jammy +dist: focal jobs: include: @@ -8,7 +8,6 @@ jobs: os: linux env: - DRAT_REPOS: "jranke" global: - USE_BSPM="true" - NOT_CRAN="true" -- cgit v1.2.1 From f21373a16d368c56234246417a32bdfb95589cf9 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Aug 2023 17:13:01 +0200 Subject: Go back to jammy on travis As it passed all tests --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 683ca4ee..3aec6949 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,6 @@ language: c sudo: required -dist: focal +dist: jammy jobs: include: -- cgit v1.2.1 From 9abab1e2d4385039b01ad3dc0d9c5966bbe94fee Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Aug 2023 17:59:42 +0200 Subject: Update static docs --- docs/404.html | 2 +- docs/articles/FOCUS_L.html | 104 +- docs/articles/index.html | 4 +- docs/articles/mkin.html | 4 +- docs/articles/prebuilt/2023_mesotrione_parent.html | 2560 ++++++++++++++++++++ .../figure-html/unnamed-chunk-14-1.png | Bin 0 -> 198820 bytes .../figure-html/unnamed-chunk-19-1.png | Bin 0 -> 195998 bytes .../figure-html/unnamed-chunk-25-1.png | Bin 0 -> 199089 bytes .../figure-html/unnamed-chunk-30-1.png | Bin 0 -> 196801 bytes .../figure-html/unnamed-chunk-8-1.png | Bin 0 -> 195237 bytes docs/authors.html | 6 +- docs/index.html | 2 +- docs/news/index.html | 8 +- docs/pkgdown.yml | 3 +- docs/reference/D24_2014.html | 2 +- docs/reference/Rplot001.png | Bin 13993 -> 88951 bytes docs/reference/Rplot002.png | Bin 13470 -> 88951 bytes docs/reference/Rplot003.png | Bin 49926 -> 88951 bytes docs/reference/dimethenamid_2018.html | 12 +- docs/reference/index.html | 2 +- docs/reference/read_spreadsheet.html | 6 +- docs/sitemap.xml | 3 + log/build.log | 2 - 23 files changed, 2645 insertions(+), 75 deletions(-) create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent.html create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png diff --git a/docs/404.html b/docs/404.html index b9ab1f25..932d8345 100644 --- a/docs/404.html +++ b/docs/404.html @@ -32,7 +32,7 @@ mkin - 1.2.4 + 1.2.5 diff --git a/docs/articles/FOCUS_L.html b/docs/articles/FOCUS_L.html index dabb16f2..7b5acf17 100644 --- a/docs/articles/FOCUS_L.html +++ b/docs/articles/FOCUS_L.html @@ -33,7 +33,7 @@ mkin - 1.2.4 + 1.2.5 @@ -135,7 +135,7 @@ to L3 Ranke

Last change 18 May 2023 -(rebuilt 2023-05-19)

+(rebuilt 2023-08-09) Source: vignettes/FOCUS_L.rmd @@ -168,18 +168,18 @@ model fit. This covers the numerical analysis given in the FOCUS report.

 m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE)
-summary(m.L1.SFO)
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:25 2023 
-## Date of summary: Fri May 19 09:20:25 2023 
+summary(m.L1.SFO)
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:39 2023 
+## Date of summary: Wed Aug  9 17:55:39 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 133 model solutions performed in 0.011 s
+## Fitted using 133 model solutions performed in 0.031 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -256,7 +256,7 @@ report.

A plot of the fit is obtained with the plot function for mkinfit objects.

-plot(m.L1.SFO, show_errmin = TRUE, main = "FOCUS L1 - SFO")
+plot(m.L1.SFO, show_errmin = TRUE, main = "FOCUS L1 - SFO")

The residual plot can be easily obtained by

@@ -268,25 +268,25 @@ objects.

## Warning in mkinfit("FOMC", FOCUS_2006_L1_mkin, quiet = TRUE): Optimisation did not converge:
 ## false convergence (8)
-plot(m.L1.FOMC, show_errmin = TRUE, main = "FOCUS L1 - FOMC")
+plot(m.L1.FOMC, show_errmin = TRUE, main = "FOCUS L1 - FOMC")

-summary(m.L1.FOMC, data = FALSE)
+summary(m.L1.FOMC, data = FALSE)
## Warning in sqrt(diag(covar)): NaNs produced
## Warning in sqrt(1/diag(V)): NaNs produced
## Warning in cov2cor(ans$covar): diag(.) had 0 or NA entries; non-finite result
 ## is doubtful
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:25 2023 
-## Date of summary: Fri May 19 09:20:25 2023 
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:39 2023 
+## Date of summary: Wed Aug  9 17:55:39 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 342 model solutions performed in 0.021 s
+## Fitted using 342 model solutions performed in 0.07 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -398,7 +398,7 @@ residual plot can be obtained simply by adding the argument
 show_residuals to the plot command.

 m.L2.SFO <- mkinfit("SFO", FOCUS_2006_L2_mkin, quiet=TRUE)
-plot(m.L2.SFO, show_residuals = TRUE, show_errmin = TRUE,
+plot(m.L2.SFO, show_residuals = TRUE, show_errmin = TRUE,
      main = "FOCUS L2 - SFO")

The \(\chi^2\) error level of 14% @@ -422,22 +422,22 @@ kinetics.

For comparison, the FOMC model is fitted as well, and the \(\chi^2\) error level is checked.

 m.L2.FOMC <- mkinfit("FOMC", FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.FOMC, show_residuals = TRUE,
+plot(m.L2.FOMC, show_residuals = TRUE,
      main = "FOCUS L2 - FOMC")

-summary(m.L2.FOMC, data = FALSE)
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:25 2023 
-## Date of summary: Fri May 19 09:20:25 2023 
+summary(m.L2.FOMC, data = FALSE)
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:40 2023 
+## Date of summary: Wed Aug  9 17:55:40 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 239 model solutions performed in 0.013 s
+## Fitted using 239 model solutions performed in 0.044 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -506,15 +506,15 @@ the data.

Fitting the four parameter DFOP model further reduces the \(\chi^2\) error level.

 m.L2.DFOP <- mkinfit("DFOP", FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
+plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
      main = "FOCUS L2 - DFOP")

-summary(m.L2.DFOP, data = FALSE)
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:25 2023 
-## Date of summary: Fri May 19 09:20:25 2023 
+summary(m.L2.DFOP, data = FALSE)
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:40 2023 
+## Date of summary: Wed Aug  9 17:55:40 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -523,7 +523,7 @@ the data.

## ## Model predictions using solution type analytical ## -## Fitted using 581 model solutions performed in 0.038 s +## Fitted using 581 model solutions performed in 0.119 s ## ## Error model: Constant variance ## @@ -611,7 +611,7 @@ only the L3 dataset prepared above.

# Only use one core here, not to offend the CRAN checks mm.L3 <- mmkin(c("SFO", "FOMC", "DFOP"), cores = 1, list("FOCUS L3" = FOCUS_2006_L3_mkin), quiet = TRUE) -plot(mm.L3)
+plot(mm.L3)

The \(\chi^2\) error level of 21% as well as the plot suggest that the SFO model does not fit very well. The @@ -627,11 +627,11 @@ as a row index and datasets as a column index.

using square brackets for indexing which will result in the use of the summary and plot functions working on mkinfit objects.

-summary(mm.L3[["DFOP", 1]])
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:26 2023 
-## Date of summary: Fri May 19 09:20:26 2023 
+summary(mm.L3[["DFOP", 1]])
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:41 2023 
+## Date of summary: Wed Aug  9 17:55:41 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -640,7 +640,7 @@ summary and plot functions working on mkinfit objects.

## ## Model predictions using solution type analytical ## -## Fitted using 376 model solutions performed in 0.021 s +## Fitted using 376 model solutions performed in 0.075 s ## ## Error model: Constant variance ## @@ -715,7 +715,7 @@ summary and plot functions working on mkinfit objects.

## 91 parent 15.0 15.18 -0.18181 ## 120 parent 12.0 10.19 1.81395
-plot(mm.L3[["DFOP", 1]], show_errmin = TRUE)
+plot(mm.L3[["DFOP", 1]], show_errmin = TRUE)

Here, a look to the model plot, the confidence intervals of the parameters and the correlation matrix suggest that the parameter @@ -746,7 +746,7 @@ below:

mm.L4 <- mmkin(c("SFO", "FOMC"), cores = 1, list("FOCUS L4" = FOCUS_2006_L4_mkin), quiet = TRUE) -plot(mm.L4)
+plot(mm.L4)

The \(\chi^2\) error level of 3.3% as well as the plot suggest that the SFO model fits very well. The error @@ -754,18 +754,18 @@ level at which the \(\chi^2\) test passes is slightly lower for the FOMC model. However, the difference appears negligible.

-summary(mm.L4[["SFO", 1]], data = FALSE)
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:26 2023 
-## Date of summary: Fri May 19 09:20:26 2023 
+summary(mm.L4[["SFO", 1]], data = FALSE)
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:42 2023 
+## Date of summary: Wed Aug  9 17:55:42 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 142 model solutions performed in 0.008 s
+## Fitted using 142 model solutions performed in 0.027 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -819,18 +819,18 @@ appears negligible.

## DT50 DT90 ## parent 106 352
-summary(mm.L4[["FOMC", 1]], data = FALSE)
-
## mkin version used for fitting:    1.2.4 
-## R version used for fitting:       4.3.0 
-## Date of fit:     Fri May 19 09:20:26 2023 
-## Date of summary: Fri May 19 09:20:26 2023 
+summary(mm.L4[["FOMC", 1]], data = FALSE)
+
## mkin version used for fitting:    1.2.5 
+## R version used for fitting:       4.3.1 
+## Date of fit:     Wed Aug  9 17:55:42 2023 
+## Date of summary: Wed Aug  9 17:55:42 2023 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 224 model solutions performed in 0.012 s
+## Fitted using 224 model solutions performed in 0.04 s
 ## 
 ## Error model: Constant variance 
 ## 
diff --git a/docs/articles/index.html b/docs/articles/index.html
index 7060b377..ac6236f1 100644
--- a/docs/articles/index.html
+++ b/docs/articles/index.html
@@ -17,7 +17,7 @@
       
       
         mkin
-        1.2.4
+        1.2.5
       
     
 
@@ -119,6 +119,8 @@
         
Testing hierarchical pathway kinetics with residue data on dimethenamid and dimethenamid-P
+
Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione
+
Calculation of time weighted average concentrations with mkin
Example evaluation of FOCUS dataset Z
diff --git a/docs/articles/mkin.html b/docs/articles/mkin.html index 6213916b..3bd0d4d3 100644 --- a/docs/articles/mkin.html +++ b/docs/articles/mkin.html @@ -33,7 +33,7 @@ mkin - 1.2.4 + 1.2.5 @@ -134,7 +134,7 @@ Ranke

Last change 18 May 2023 -(rebuilt 2023-05-19)

+(rebuilt 2023-08-09) Source: vignettes/mkin.rmd diff --git a/docs/articles/prebuilt/2023_mesotrione_parent.html b/docs/articles/prebuilt/2023_mesotrione_parent.html new file mode 100644 index 00000000..b233fc3c --- /dev/null +++ b/docs/articles/prebuilt/2023_mesotrione_parent.html @@ -0,0 +1,2560 @@ + + + + + + + +Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione • mkin + + + + + + + + + + + + +
+
+ + + + +
+
+ + + + +
+

Introduction +

+

The purpose of this document is to test demonstrate how nonlinear +hierarchical models (NLHM) based on the parent degradation models SFO, +FOMC, DFOP and HS can be fitted with the mkin package, also considering +the influence of covariates like soil pH on different degradation +parameters. Because in some other case studies, the SFORB +parameterisation of biexponential decline has shown some advantages over +the DFOP parameterisation, SFORB was included in the list of tested +models as well.

+

The mkin package is used in version 1.2.5, which is contains the +functions that were used for the evaluations. The saemix +package is used as a backend for fitting the NLHM, but is also loaded to +make the convergence plot function available.

+

This document is processed with the knitr package, which +also provides the kable function that is used to improve +the display of tabular data in R markdown documents. For parallel +processing, the parallel package is used.

+
+library(mkin)
+library(knitr)
+library(saemix)
+library(parallel)
+n_cores <- detectCores()
+if (Sys.info()["sysname"] == "Windows") {
+  cl <- makePSOCKcluster(n_cores)
+} else {
+  cl <- makeForkCluster(n_cores)
+}
+
+

Test data +

+
+data_file <- system.file(
+  "testdata", "mesotrione_soil_efsa_2016.xlsx", package = "mkin")
+meso_ds <- read_spreadsheet(data_file, parent_only = TRUE)
+

The following tables show the covariate data and the 18 datasets that +were read in from the spreadsheet file.

+
+pH <- attr(meso_ds, "covariates")
+kable(pH, caption = "Covariate data")
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Covariate data
pH
Richmond6.2
Richmond 26.2
ERTC6.4
Toulouse7.7
Picket Piece7.1
7215.6
7225.7
7235.4
7244.8
7255.8
7275.1
7285.9
7295.6
7305.3
7316.1
7325.0
7415.7
7427.2
+
+for (ds_name in names(meso_ds)) {
+  print(
+    kable(mkin_long_to_wide(meso_ds[[ds_name]]),
+      caption = paste("Dataset", ds_name),
+      booktabs = TRUE, row.names = FALSE))
+}
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Richmond
timemeso
0.00000091.00
1.17905086.70
3.53714973.60
7.07429961.50
10.61144855.70
15.32764747.70
17.68574739.50
24.76004629.80
35.37149419.60
68.3848895.67
0.00000097.90
1.17905096.40
3.53714989.10
7.07429974.40
10.61144857.40
15.32764746.30
18.86479735.50
27.11814627.20
35.37149419.10
74.2801386.50
108.4725823.40
142.6650272.20
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Richmond 2
timemeso
0.00000096.0
2.42200482.4
5.65134371.2
8.07334853.1
11.30268748.5
16.95403033.4
22.60537324.2
45.21074611.9
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset ERTC
timemeso
0.00000099.9
2.75519380.0
6.42878242.1
9.18397550.1
12.85756528.4
19.28634739.8
25.71513029.9
51.4302592.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Toulouse
timemeso
0.00000096.8
2.89798363.3
6.76196022.3
9.65994216.6
13.52391916.1
20.28587917.2
27.0478381.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Picket Piece
timemeso
0.000000102.0
2.84119573.7
6.62945435.5
9.47064931.8
13.25890918.0
19.8883643.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 721
timemeso
0.0000086.4
11.2436661.4
22.4873349.8
33.7309941.0
44.9746635.1
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 722
timemeso
0.0000090.3
11.2436652.1
22.4873337.4
33.7309921.2
44.9746614.3
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 723
timemeso
0.0000089.3
11.2436670.8
22.4873351.1
33.7309942.7
44.9746626.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 724
timemeso
0.00000089.4
9.00820865.2
18.01641555.8
27.02462346.0
36.03283141.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 725
timemeso
0.0000089.0
10.9905835.4
21.9811618.6
32.9717411.6
43.962327.6
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 727
timemeso
0.0000091.3
10.9610463.2
21.9220951.1
32.8831342.0
43.8441740.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 728
timemeso
0.0000091.8
11.2436643.6
22.4873322.0
33.7309915.9
44.974668.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 729
timemeso
0.0000091.6
11.2436660.5
22.4873343.5
33.7309928.4
44.9746620.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 730
timemeso
0.0000092.7
11.0744658.9
22.1489344.0
33.2233946.0
44.2978529.3
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 731
timemeso
0.0000092.1
11.2436664.4
22.4873345.3
33.7309933.6
44.9746623.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 732
timemeso
0.0000090.3
11.2436658.2
22.4873340.1
33.7309933.1
44.9746625.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 741
timemeso
0.0000090.3
10.8471268.7
21.6942458.0
32.5413652.2
43.3884848.0
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 742
timemeso
0.0000092.0
11.2436660.9
22.4873336.2
33.7309918.3
44.974668.7
+
+
+
+

Separate evaluations +

+

In order to obtain suitable starting parameters for the NLHM fits, +separate fits of the five models to the data for each soil are generated +using the mmkin function from the mkin package. In a first +step, constant variance is assumed. Convergence is checked with the +status function.

+
+deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB", "HS")
+f_sep_const <- mmkin(
+  deg_mods,
+  meso_ds,
+  error_model = "const",
+  cluster = cl,
+  quiet = TRUE)
+
+status(f_sep_const[, 1:5]) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
RichmondRichmond 2ERTCToulousePicket Piece
SFOOKOKOKOKOK
FOMCOKOKOKOKC
DFOPOKOKOKOKOK
SFORBOKOKOKOKOK
HSOKOKCOKOK
+
+status(f_sep_const[, 6:18]) |> kable()
+ ++++++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
721722723724725727728729730731732741742
SFOOKOKOKOKOKOKOKOKOKOKOKOKOK
FOMCOKOKCOKOKOKOKOKOKOKOKOKOK
DFOPOKOKOKOKOKOKOKOKOKOKOKOKOK
SFORBOKOKOKOKOKOKOKCOKOKOKOKOK
HSOKOKOKOKOKOKOKOKOKOKOKOKOK
+

In the tables above, OK indicates convergence and C indicates failure +to converge. Most separate fits with constant variance converged, with +the exception of two FOMC fits, one SFORB fit and one HS fit.

+
+f_sep_tc <- update(f_sep_const, error_model = "tc")
+
+status(f_sep_tc[, 1:5]) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
RichmondRichmond 2ERTCToulousePicket Piece
SFOOKOKOKOKOK
FOMCOKOKOKOKOK
DFOPCOKOKOKOK
SFORBOKOKOKOKOK
HSOKOKCOKOK
+
+status(f_sep_tc[, 6:18]) |> kable()
+ ++++++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
721722723724725727728729730731732741742
SFOOKOKOKOKOKOKOKOKOKOKOKOKOK
FOMCOKOKCOKCCOKCOKCOKCOK
DFOPCOKOKOKCOKOKOKOKCOKCOK
SFORBCOKOKOKCOKOKCOKOKOKCOK
HSOKOKOKOKOKOKOKOKOKCOKOKOK
+

With the two-component error model, the set of fits that did not +converge is larger, with convergence problems appearing for a number of +non-SFO fits.

+
+
+

Hierarchical model fits without covariate effect +

+

The following code fits hierarchical kinetic models for the ten +combinations of the five different degradation models with the two +different error models in parallel.

+
+f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cluster = cl)
+status(f_saem_1) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOOKOK
FOMCOKOK
DFOPOKOK
SFORBOKOK
HSOKOK
+

All fits terminate without errors (status OK).

+
+anova(f_saem_1) |> kable(digits = 1)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
nparAICBICLik
SFO const5800.0804.5-395.0
SFO tc6801.9807.2-394.9
FOMC const7787.4793.6-386.7
FOMC tc8788.9796.1-386.5
DFOP const9787.6795.6-384.8
SFORB const9787.4795.4-384.7
HS const9781.9789.9-382.0
DFOP tc10787.4796.3-383.7
SFORB tc10795.8804.7-387.9
HS tc10783.7792.7-381.9
+

The model comparisons show that the fits with constant variance are +consistently preferable to the corresponding fits with two-component +error for these data. This is confirmed by the fact that the parameter +b.1 (the relative standard deviation in the fits obtained +with the saemix package), is ill-defined in all fits.

+
+illparms(f_saem_1) |> kable()
+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOsd(meso_0)sd(meso_0), b.1
FOMCsd(meso_0), sd(log_beta)sd(meso_0), sd(log_beta), b.1
DFOPsd(meso_0), sd(log_k1)sd(meso_0), sd(g_qlogis), b.1
SFORBsd(meso_free_0), sd(log_k_meso_free_bound)sd(meso_free_0), sd(log_k_meso_free_bound), b.1
HSsd(meso_0)sd(meso_0), b.1
+

For obtaining fits with only well-defined random effects, we update +the set of fits, excluding random effects that were ill-defined +according to the illparms function.

+
+f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+status(f_saem_2) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOOKOK
FOMCOKOK
DFOPOKOK
SFORBOKOK
HSOKOK
+

The updated fits terminate without errors.

+
+illparms(f_saem_2) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOb.1
FOMCb.1
DFOPb.1
SFORBb.1
HS
+

No ill-defined errors remain in the fits with constant variance.

+
+
+

Hierarchical model fits with covariate effect +

+

In the following sections, hierarchical fits including a model for +the influence of pH on selected degradation parameters are shown for all +parent models. Constant variance is selected as the error model based on +the fits without covariate effects. Random effects that were ill-defined +in the fits without pH influence are excluded. A potential influence of +the soil pH is only included for parameters with a well-defined random +effect, because experience has shown that only for such parameters a +significant pH effect could be found.

+
+

SFO +

+
+sfo_pH <- saem(f_sep_const["SFO", ], no_random_effect = "meso_0", covariates = pH,
+  covariate_models = list(log_k_meso ~ pH))
+
+summary(sfo_pH)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_091.3589.2793.43
log_k_meso-6.66-7.97-5.35
beta_pH(log_k_meso)0.590.370.81
a.15.484.716.24
SD.log_k_meso0.350.230.47
+

The parameter showing the pH influence in the above table is +beta_pH(log_k_meso). Its confidence interval does not +include zero, indicating that the influence of soil pH on the log of the +degradation rate constant is significantly greater than zero.

+
+anova(f_saem_2[["SFO", "const"]], sfo_pH, test = TRUE)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+                           npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
+f_saem_2[["SFO", "const"]]    4 797.56 801.12 -394.78                         
+sfo_pH                        5 783.09 787.54 -386.54 16.473  1  4.934e-05 ***
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+

The comparison with the SFO fit without covariate effect confirms +that considering the soil pH improves the model, both by comparison of +AIC and BIC and by the likelihood ratio test.

+
+plot(sfo_pH)
+

+

Endpoints for a model with covariates are by default calculated for +the median of the covariate values. This quantile can be adapted, or a +specific covariate value can be given as shown below.

+
+endpoints(sfo_pH)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90
+meso 18.52069 61.52441
+
+endpoints(sfo_pH, covariate_quantile = 0.9)
+
$covariates
+      pH
+90% 7.13
+
+$distimes
+         DT50     DT90
+meso 8.237019 27.36278
+
+endpoints(sfo_pH, covariates = c(pH = 7.0))
+
$covariates
+     pH
+User  7
+
+$distimes
+        DT50    DT90
+meso 8.89035 29.5331
+
+
+

FOMC +

+
+fomc_pH <- saem(f_sep_const["FOMC", ], no_random_effect = "meso_0", covariates = pH,
+  covariate_models = list(log_alpha ~ pH))
+
+summary(fomc_pH)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_092.8490.7594.93
log_alpha-2.21-3.49-0.92
beta_pH(log_alpha)0.580.370.79
log_beta4.213.444.99
a.15.034.325.73
SD.log_alpha0.00-23.7723.78
SD.log_beta0.370.010.74
+

As in the case of SFO, the confidence interval of the slope parameter +(here beta_pH(log_alpha)) quantifying the influence of soil +pH does not include zero, and the model comparison clearly indicates +that the model with covariate influence is preferable. However, the +random effect for alpha is not well-defined any more after +inclusion of the covariate effect (the confidence interval of +SD.log_alpha includes zero).

+
+illparms(fomc_pH)
+
[1] "sd(log_alpha)"
+

Therefore, the model is updated without this random effect, and no +ill-defined parameters remain.

+
+fomc_pH_2 <- update(fomc_pH, no_random_effect = c("meso_0", "log_alpha"))
+illparms(fomc_pH_2)
+
+anova(f_saem_2[["FOMC", "const"]], fomc_pH, fomc_pH_2, test = TRUE)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+                            npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
+f_saem_2[["FOMC", "const"]]    5 783.25 787.71 -386.63                         
+fomc_pH_2                      6 767.49 772.83 -377.75 17.762  1  2.503e-05 ***
+fomc_pH                        7 770.07 776.30 -378.04  0.000  1          1    
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+

Model comparison indicates that including pH dependence significantly +improves the fit, and that the reduced model with covariate influence +results in the most preferable FOMC fit.

+
+summary(fomc_pH_2)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.0590.9895.13
log_alpha-2.91-4.18-1.63
beta_pH(log_alpha)0.660.440.87
log_beta3.953.294.62
a.14.984.285.68
SD.log_beta0.400.260.54
+
+plot(fomc_pH_2)
+

+
+endpoints(fomc_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90 DT50back
+meso 17.30248 82.91343 24.95943
+
+endpoints(fomc_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back
+meso 6.986239 27.02927 8.136621
+
+
+

DFOP +

+

In the DFOP fits without covariate effects, random effects for two +degradation parameters (k2 and g) were +identifiable.

+
+summary(f_saem_2[["DFOP", "const"]])$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.6191.5895.63
log_k1-1.53-2.27-0.79
log_k2-3.42-3.73-3.11
g_qlogis-1.67-2.57-0.77
a.14.744.025.45
SD.log_k20.600.380.81
SD.g_qlogis0.940.331.54
+

A fit with pH dependent degradation parameters was obtained by +excluding the same random effects as in the refined DFOP fit without +covariate influence, and including covariate models for the two +identifiable parameters k2 and g.

+
+dfop_pH <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
+  covariates = pH,
+  covariate_models = list(log_k2 ~ pH, g_qlogis ~ pH))
+

The corresponding parameters for the influence of soil pH are +beta_pH(log_k2) for the influence of soil pH on +k2, and beta_pH(g_qlogis) for its influence on +g.

+
+summary(dfop_pH)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_092.8490.8594.84
log_k1-2.82-3.09-2.54
log_k2-11.48-15.32-7.64
beta_pH(log_k2)1.310.691.92
g_qlogis3.130.475.80
beta_pH(g_qlogis)-0.57-1.04-0.09
a.14.964.265.65
SD.log_k20.760.471.05
SD.g_qlogis0.01-9.969.97
+
+illparms(dfop_pH)
+
[1] "sd(g_qlogis)"
+

Confidence intervals for neither of them include zero, indicating a +significant difference from zero. However, the random effect for +g is now ill-defined. The fit is updated without this +ill-defined random effect.

+
+dfop_pH_2 <- update(dfop_pH,
+  no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
+illparms(dfop_pH_2)
+
[1] "beta_pH(g_qlogis)"
+

Now, the slope parameter for the pH effect on g is +ill-defined. Therefore, another attempt is made without the +corresponding covariate model.

+
+dfop_pH_3 <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
+  covariates = pH,
+  covariate_models = list(log_k2 ~ pH))
+illparms(dfop_pH_3)
+
[1] "sd(g_qlogis)"
+

As the random effect for g is again ill-defined, the fit +is repeated without it.

+
+dfop_pH_4 <- update(dfop_pH_3, no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
+illparms(dfop_pH_4)
+

While no ill-defined parameters remain, model comparison suggests +that the previous model dfop_pH_2 with two pH dependent +parameters is preferable, based on information criteria as well as based +on the likelihood ratio test.

+
+anova(f_saem_2[["DFOP", "const"]], dfop_pH, dfop_pH_2, dfop_pH_3, dfop_pH_4)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+                            npar    AIC    BIC     Lik
+f_saem_2[["DFOP", "const"]]    7 782.94 789.18 -384.47
+dfop_pH_4                      7 767.35 773.58 -376.68
+dfop_pH_2                      8 765.14 772.26 -374.57
+dfop_pH_3                      8 769.00 776.12 -376.50
+dfop_pH                        9 769.10 777.11 -375.55
+
+anova(dfop_pH_2, dfop_pH_4, test = TRUE)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+          npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)  
+dfop_pH_4    7 767.35 773.58 -376.68                       
+dfop_pH_2    8 765.14 772.26 -374.57 4.2153  1    0.04006 *
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+

When focussing on parameter identifiability using the test if the +confidence interval includes zero, dfop_pH_4 would still be +the preferred model. However, it should be kept in mind that parameter +confidence intervals are constructed using a simple linearisation of the +likelihood. As the confidence interval of the random effect for +g only marginally includes zero, it is suggested that this +is acceptable, and that dfop_pH_2 can be considered the +most preferable model.

+
+plot(dfop_pH_2)
+

+
+endpoints(dfop_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 18.36876 73.51841 22.13125 4.191901 23.98672
+
+endpoints(dfop_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 8.346428 28.34437 8.532507 4.191901 8.753618
+
+
+

SFORB +

+
+sforb_pH <- saem(f_sep_const["SFORB", ], no_random_effect = c("meso_free_0", "log_k_meso_free_bound"),
+  covariates = pH,
+  covariate_models = list(log_k_meso_free ~ pH, log_k_meso_bound_free ~ pH))
+
+summary(sforb_pH)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_free_093.4291.3295.52
log_k_meso_free-5.37-6.94-3.81
beta_pH(log_k_meso_free)0.420.180.67
log_k_meso_free_bound-3.49-4.92-2.05
log_k_meso_bound_free-9.98-19.22-0.74
beta_pH(log_k_meso_bound_free)1.23-0.212.67
a.14.904.185.63
SD.log_k_meso_free0.350.230.47
SD.log_k_meso_bound_free0.13-1.952.20
+

The confidence interval of +beta_pH(log_k_meso_bound_free) includes zero, indicating +that the influence of soil pH on k_meso_bound_free cannot +reliably be quantified. Also, the confidence interval for the random +effect on this parameter (SD.log_k_meso_bound_free) +includes zero.

+

Using the illparms function, these ill-defined +parameters can be found more conveniently.

+
+illparms(sforb_pH)
+
[1] "sd(log_k_meso_bound_free)"      "beta_pH(log_k_meso_bound_free)"
+

To remove the ill-defined parameters, a second variant of the SFORB +model with pH influence is fitted. No ill-defined parameters remain.

+
+sforb_pH_2 <- update(sforb_pH,
+  no_random_effect = c("meso_free_0", "log_k_meso_free_bound", "log_k_meso_bound_free"),
+  covariate_models = list(log_k_meso_free ~ pH))
+illparms(sforb_pH_2)
+

The model comparison of the SFORB fits includes the refined model +without covariate effect, and both versions of the SFORB fit with +covariate effect.

+
+anova(f_saem_2[["SFORB", "const"]], sforb_pH, sforb_pH_2, test = TRUE)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+                             npar    AIC    BIC     Lik   Chisq Df Pr(>Chisq)  
+f_saem_2[["SFORB", "const"]]    7 783.40 789.63 -384.70                        
+sforb_pH_2                      7 770.94 777.17 -378.47 12.4616  0             
+sforb_pH                        9 768.81 776.83 -375.41  6.1258  2    0.04675 *
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+

The first model including pH influence is preferable based on +information criteria and the likelihood ratio test. However, as it is +not fully identifiable, the second model is selected.

+
+summary(sforb_pH_2)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_free_093.3291.1695.48
log_k_meso_free-6.15-7.43-4.86
beta_pH(log_k_meso_free)0.540.330.75
log_k_meso_free_bound-3.80-5.20-2.40
log_k_meso_bound_free-2.95-4.26-1.64
a.15.084.385.79
SD.log_k_meso_free0.330.220.45
+
+plot(sforb_pH_2)
+

+
+endpoints(sforb_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$ff
+meso_free 
+        1 
+
+$SFORB
+   meso_b1    meso_b2     meso_g 
+0.09735824 0.02631699 0.31602120 
+
+$distimes
+         DT50     DT90 DT50back DT50_meso_b1 DT50_meso_b2
+meso 16.86549 73.15824 22.02282     7.119554     26.33839
+
+endpoints(sforb_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$ff
+meso_free 
+        1 
+
+$SFORB
+   meso_b1    meso_b2     meso_g 
+0.13315233 0.03795988 0.61186191 
+
+$distimes
+         DT50     DT90 DT50back DT50_meso_b1 DT50_meso_b2
+meso 7.932495 36.93311 11.11797     5.205671        18.26
+
+
+

HS +

+
+hs_pH <- saem(f_sep_const["HS", ], no_random_effect = c("meso_0"),
+  covariates = pH,
+  covariate_models = list(log_k1 ~ pH, log_k2 ~ pH, log_tb ~ pH))
+
+summary(hs_pH)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.3391.4795.19
log_k1-5.81-7.27-4.36
beta_pH(log_k1)0.470.230.72
log_k2-6.80-8.76-4.83
beta_pH(log_k2)0.540.210.87
log_tb3.251.255.25
beta_pH(log_tb)-0.10-0.430.23
a.14.493.785.21
SD.log_k10.370.240.51
SD.log_k20.290.100.48
SD.log_tb0.25-0.070.57
+
+illparms(hs_pH)
+
[1] "sd(log_tb)"      "beta_pH(log_tb)"
+

According to the output of the illparms function, the +random effect on the break time tb cannot reliably be +quantified, neither can the influence of soil pH on tb. The +fit is repeated without the corresponding covariate model, and no +ill-defined parameters remain.

+
+hs_pH_2 <- update(hs_pH, covariate_models = list(log_k1 ~ pH, log_k2 ~ pH))
+illparms(hs_pH_2)
+

Model comparison confirms that this model is preferable to the fit +without covariate influence, and also to the first version with +covariate influence.

+
+anova(f_saem_2[["HS", "const"]], hs_pH, hs_pH_2, test = TRUE)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+                          npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
+f_saem_2[["HS", "const"]]    8 780.08 787.20 -382.04                         
+hs_pH_2                     10 766.47 775.37 -373.23 17.606  2  0.0001503 ***
+hs_pH                       11 769.80 779.59 -373.90  0.000  1  1.0000000    
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
+summary(hs_pH_2)$confint_trans |> kable(digits = 2)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.3391.5095.15
log_k1-5.68-7.09-4.27
beta_pH(log_k1)0.460.220.69
log_k2-6.61-8.34-4.88
beta_pH(log_k2)0.500.210.79
log_tb2.702.333.08
a.14.453.745.16
SD.log_k10.360.220.49
SD.log_k20.230.020.43
SD.log_tb0.550.250.85
+
+plot(hs_pH_2)
+

+
+endpoints(hs_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 14.68725 82.45287 24.82079 14.68725 29.29299
+
+endpoints(hs_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 8.298536 38.85371 11.69613 8.298536 15.71561
+
+
+

Comparison across parent models +

+

After model reduction for all models with pH influence, they are +compared with each other.

+
+anova(sfo_pH, fomc_pH_2, dfop_pH_2, dfop_pH_4, sforb_pH_2, hs_pH_2)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+           npar    AIC    BIC     Lik
+sfo_pH        5 783.09 787.54 -386.54
+fomc_pH_2     6 767.49 772.83 -377.75
+dfop_pH_4     7 767.35 773.58 -376.68
+sforb_pH_2    7 770.94 777.17 -378.47
+dfop_pH_2     8 765.14 772.26 -374.57
+hs_pH_2      10 766.47 775.37 -373.23
+

The DFOP model with pH influence on k2 and +g and a random effect only on k2 is finally +selected as the best fit.

+

The endpoints resulting from this model are listed below. Please +refer to the Appendix for a detailed listing.

+
+endpoints(dfop_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 18.36876 73.51841 22.13125 4.191901 23.98672
+
+endpoints(dfop_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 8.346428 28.34437 8.532507 4.191901 8.753618
+
+
+
+

Conclusions +

+

These evaluations demonstrate that covariate effects can be included +for all types of parent degradation models. These models can then be +further refined to make them fully identifiable.

+
+
+

Appendix +

+
+

Hierarchical fit listings +

+
+

Fits without covariate effects +

+ +
+
+

Fits with covariate effects +

+ +
+
+
+

Session info +

+
R version 4.3.1 (2023-06-16)
+Platform: x86_64-pc-linux-gnu (64-bit)
+Running under: Debian GNU/Linux 12 (bookworm)
+
+Matrix products: default
+BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
+LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0
+
+locale:
+ [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
+ [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
+ [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
+ [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
+ [9] LC_ADDRESS=C               LC_TELEPHONE=C            
+[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
+
+time zone: Europe/Berlin
+tzcode source: system (glibc)
+
+attached base packages:
+[1] parallel  stats     graphics  grDevices utils     datasets  methods  
+[8] base     
+
+other attached packages:
+[1] saemix_3.2 npde_3.3   knitr_1.43 mkin_1.2.5
+
+loaded via a namespace (and not attached):
+ [1] sass_0.4.6        utf8_1.2.3        generics_0.1.3    stringi_1.7.12   
+ [5] lattice_0.20-45   digest_0.6.31     magrittr_2.0.3    evaluate_0.21    
+ [9] grid_4.3.1        fastmap_1.1.1     cellranger_1.1.0  rprojroot_2.0.3  
+[13] jsonlite_1.8.5    mclust_6.0.0      gridExtra_2.3     purrr_1.0.1      
+[17] fansi_1.0.4       scales_1.2.1      codetools_0.2-19  textshaping_0.3.6
+[21] jquerylib_0.1.4   cli_3.6.1         rlang_1.1.1       munsell_0.5.0    
+[25] cachem_1.0.8      yaml_2.3.7        tools_4.3.1       memoise_2.0.1    
+[29] dplyr_1.1.2       colorspace_2.1-0  ggplot2_3.4.2     vctrs_0.6.2      
+[33] R6_2.5.1          zoo_1.8-12        lifecycle_1.0.3   stringr_1.5.0    
+[37] fs_1.6.2          ragg_1.2.5        pkgconfig_2.0.3   desc_1.4.2       
+[41] pkgdown_2.0.7     bslib_0.4.2       pillar_1.9.0      gtable_0.3.3     
+[45] glue_1.6.2        systemfonts_1.0.4 highr_0.10        xfun_0.39        
+[49] tibble_3.2.1      lmtest_0.9-40     tidyselect_1.2.0  htmltools_0.5.5  
+[53] nlme_3.1-162      rmarkdown_2.22    compiler_4.3.1    readxl_1.4.2     
+
+
+

Hardware info +

+
CPU model: Intel(R) Core(TM) i7-4710MQ CPU @ 2.50GHz
+
MemTotal:       12165632 kB
+
+
+
+ + + +
+ + + +
+ +
+

+

Site built with pkgdown 2.0.7.

+
+ +
+
+ + + + + + + + diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png new file mode 100644 index 00000000..863a48bd Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png new file mode 100644 index 00000000..256b2b68 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png new file mode 100644 index 00000000..59011020 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png new file mode 100644 index 00000000..f427bc39 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png new file mode 100644 index 00000000..7c3b460b Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/docs/authors.html b/docs/authors.html index 313c86e9..8b91ee77 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -17,7 +17,7 @@ mkin - 1.2.4 + 1.2.5 @@ -132,13 +132,13 @@

Ranke J (2023). mkin: Kinetic Evaluation of Chemical Degradation Data. -R package version 1.2.4, https://pkgdown.jrwb.de/mkin/. +R package version 1.2.5, https://pkgdown.jrwb.de/mkin/.

@Manual{,
   title = {mkin: Kinetic Evaluation of Chemical Degradation Data},
   author = {Johannes Ranke},
   year = {2023},
-  note = {R package version 1.2.4},
+  note = {R package version 1.2.5},
   url = {https://pkgdown.jrwb.de/mkin/},
 }
diff --git a/docs/index.html b/docs/index.html index 06cdfdfb..79716087 100644 --- a/docs/index.html +++ b/docs/index.html @@ -44,7 +44,7 @@ mkin - 1.2.4 + 1.2.5 diff --git a/docs/news/index.html b/docs/news/index.html index 16fce355..7f7c0a6d 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -17,7 +17,7 @@ mkin - 1.2.4 + 1.2.5 @@ -104,6 +104,12 @@ Source: NEWS.md +
+ +
  • ‘vignettes/mesotrione_parent_2023.rnw’: Prebuilt vignette showing how covariate modelling can be done for all relevant parent degradation models.

  • +
  • ‘inst/testdata/mesotrione_soil_efsa_2016}.xlsx’: Another example spreadsheets for use with ‘read_spreadsheet()’, featuring pH dependent degradation

  • +
  • R/illparms.R: Fix the detection of ill-defined slope or error model parameters for the case that the estimate is negative

  • +
  • R/endpoints.R: Fix the calculation of endpoints for user specified covariate values
  • diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index fdcd875b..65d77082 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -8,6 +8,7 @@ articles: 2022_cyan_pathway: prebuilt/2022_cyan_pathway.html 2022_dmta_parent: prebuilt/2022_dmta_parent.html 2022_dmta_pathway: prebuilt/2022_dmta_pathway.html + 2023_mesotrione_parent: prebuilt/2023_mesotrione_parent.html twa: twa.html FOCUS_Z: web_only/FOCUS_Z.html NAFTA_examples: web_only/NAFTA_examples.html @@ -16,7 +17,7 @@ articles: dimethenamid_2018: web_only/dimethenamid_2018.html multistart: web_only/multistart.html saem_benchmarks: web_only/saem_benchmarks.html -last_built: 2023-05-19T09:04Z +last_built: 2023-08-09T15:42Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/reference/D24_2014.html b/docs/reference/D24_2014.html index fcc3cec3..1e35e864 100644 --- a/docs/reference/D24_2014.html +++ b/docs/reference/D24_2014.html @@ -22,7 +22,7 @@ constrained by data protection regulations.">