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