Introduction

+ + +

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, with serial formation of two or more metabolites can @@ -377,7 +103,7 @@ be fitted with the mkin package.

173340 (Application of nonlinear hierarchical models to the kinetic evaluation of chemical degradation data) of the German Environment Agency carried out in 2022 and 2023.

The mkin package is used in version 1.2.9 which is currently under +

The mkin package is used in version 1.2.10 which is currently under development. The newly introduced functionality that is used here is a simplification of excluding random effects for a set of fits based on a related set of fits with a reduced model, and the documentation of the @@ -389,47 +115,49 @@ but is also loaded to make the convergence plot function available.

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()
-
-# We need to start a new cluster after defining a compiled model that is
-# saved as a DLL to the user directory, therefore we define a function
-# This is used again after defining the pathway model
-start_cluster <- function(n_cores) {
-  if (Sys.info()["sysname"] == "Windows") {
-    ret <- makePSOCKcluster(n_cores)
-  } else {
-    ret <- makeForkCluster(n_cores)
-  }
-  return(ret)
-}
-cl <- start_cluster(n_cores)

Test data

+library(mkin)
+library(knitr)
+library(saemix)
+library(parallel)
+n_cores <- detectCores()
+
+# We need to start a new cluster after defining a compiled model that is
+# saved as a DLL to the user directory, therefore we define a function
+# This is used again after defining the pathway model
+start_cluster <- function(n_cores) {
+  if (Sys.info()["sysname"] == "Windows") {
+    ret <- makePSOCKcluster(n_cores)
+  } else {
+    ret <- makeForkCluster(n_cores)
+  }
+  return(ret)
+}
+cl <- start_cluster(n_cores)

Test data +

The example data are taken from the final addendum to the DAR from 2014 and are distributed with the mkin package. Residue data and time step normalisation factors are read in using the function read_spreadsheet from the mkin package. This function also performs the time step normalisation.

data_file <- system.file(
-  "testdata", "cyantraniliprole_soil_efsa_2014.xlsx",
-  package = "mkin")
-cyan_ds <- read_spreadsheet(data_file, parent_only = FALSE)

+data_file <- system.file(
+  "testdata", "cyantraniliprole_soil_efsa_2014.xlsx",
+  package = "mkin")
+cyan_ds <- read_spreadsheet(data_file, parent_only = FALSE)

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

pH <- attr(cyan_ds, "covariates")
-kable(pH, caption = "Covariate data")

- +

+pH <- attr(cyan_ds, "covariates")
+kable(pH, caption = "Covariate data")

- - + - - + @@ -453,25 +181,24 @@ kable(pH, caption = "Covariate data")

Covariate data
	pH
Nambsheim

for (ds_name in names(cyan_ds)) {
-  print(
-    kable(mkin_long_to_wide(cyan_ds[[ds_name]]),
-      caption = paste("Dataset", ds_name),
-      booktabs = TRUE, row.names = FALSE))
-    cat("\n\\clearpage\n")
-}

- +

+for (ds_name in names(cyan_ds)) {
+  print(
+    kable(mkin_long_to_wide(cyan_ds[[ds_name]]),
+      caption = paste("Dataset", ds_name),
+      booktabs = TRUE, row.names = FALSE))
+    cat("\n\\clearpage\n")
+}

- - + - - + @@ -683,17 +410,15 @@ kable(pH, caption = "Covariate data")

Dataset Nambsheim
time	cyan	JCZ38	J9C38	JSE76	J9Z38
0.000000

- +

- - + - - + @@ -879,17 +604,15 @@ kable(pH, caption = "Covariate data")

Dataset Tama
time	cyan	JCZ38	J9Z38	JSE76
0.000000

- +

- - + - - + @@ -1145,17 +868,15 @@ kable(pH, caption = "Covariate data")

Dataset Gross-Umstadt
time	cyan	JCZ38	J9Z38	JSE76
0.0000000

- +

- - + - - + @@ -1285,17 +1006,15 @@ kable(pH, caption = "Covariate data")

Dataset Sassafras
time	cyan	JCZ38	J9Z38	JSE76
0.000000

- +

- - + - - + @@ -1427,25 +1146,25 @@ kable(pH, caption = "Covariate data")

Dataset Lleida
time	cyan	JCZ38	J9Z38	JSE76
0.000000

Parent only evaluations

Parent only evaluations +

As the pathway fits have very long run times, evaluations of the parent data are performed first, in order to determine for each hierarchical parent degradation model which random effects on the degradation model parameters are ill-defined.

cyan_sep_const <- mmkin(c("SFO", "FOMC", "DFOP", "SFORB", "HS"),
-  cyan_ds, quiet = TRUE, cores = n_cores)
-cyan_sep_tc <- update(cyan_sep_const, error_model = "tc")
-cyan_saem_full <- mhmkin(list(cyan_sep_const, cyan_sep_tc))
-status(cyan_saem_full) |> kable()

- - - +

+cyan_sep_const <- mmkin(c("SFO", "FOMC", "DFOP", "SFORB", "HS"),
+  cyan_ds, quiet = TRUE, cores = n_cores)
+cyan_sep_tc <- update(cyan_sep_const, error_model = "tc")
+cyan_saem_full <- mhmkin(list(cyan_sep_const, cyan_sep_tc))
+status(cyan_saem_full) |> kable()

+ - - + @@ -1475,15 +1194,14 @@ status(cyan_saem_full) |> kable()

	const	tc
SFO

All fits converged successfully.

illparms(cyan_saem_full) |> kable()

- - - +

+illparms(cyan_saem_full) |> kable()

+ - - + @@ -1516,17 +1234,16 @@ status(cyan_saem_full) |> kable() of the parent compound is ill-defined. For the biexponential models DFOP and SFORB, the random effect of one additional parameter is ill-defined when the two-component error model is used.

anova(cyan_saem_full) |> kable(digits = 1)

	const	tc
SFO

- - +

+anova(cyan_saem_full) |> kable(digits = 1)

+ - - + @@ -1604,71 +1321,74 @@ when the two-component error model is used.

two-component error model is preferable for all parent models with the exception of DFOP. The lowest AIC and BIC values are are obtained with the FOMC model, followed by SFORB and DFOP.

stopCluster(cl)

+stopCluster(cl)

Pathway fits

Evaluations with pathway established previously

Pathway fits +

Evaluations with pathway established previously +

To test the technical feasibility of coupling the relevant parent degradation models with different transformation pathway models, a list of mkinmod models is set up below. As in the EU evaluation, parallel formation of metabolites JCZ38 and J9Z38 and secondary formation of metabolite JSE76 from JCZ38 is used.

if (!dir.exists("cyan_dlls")) dir.create("cyan_dlls")
-cyan_path_1 <- list(
-  sfo_path_1 = mkinmod(
-    cyan = mkinsub("SFO", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO"), quiet = TRUE,
-    name = "sfo_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
-  fomc_path_1 = mkinmod(
-    cyan = mkinsub("FOMC", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO"), quiet = TRUE,
-    name = "fomc_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
-  dfop_path_1 = mkinmod(
-    cyan = mkinsub("DFOP", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO"), quiet = TRUE,
-    name = "dfop_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
-  sforb_path_1 = mkinmod(
-    cyan = mkinsub("SFORB", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO"), quiet = TRUE,
-    name = "sforb_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
-  hs_path_1 = mkinmod(
-    cyan = mkinsub("HS", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO"), quiet = TRUE,
-    name = "hs_path_1", dll_dir = "cyan_dlls", overwrite = TRUE)
-)
-cl_path_1 <- start_cluster(n_cores)

+if (!dir.exists("cyan_dlls")) dir.create("cyan_dlls")
+cyan_path_1 <- list(
+  sfo_path_1 = mkinmod(
+    cyan = mkinsub("SFO", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO"), quiet = TRUE,
+    name = "sfo_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
+  fomc_path_1 = mkinmod(
+    cyan = mkinsub("FOMC", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO"), quiet = TRUE,
+    name = "fomc_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
+  dfop_path_1 = mkinmod(
+    cyan = mkinsub("DFOP", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO"), quiet = TRUE,
+    name = "dfop_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
+  sforb_path_1 = mkinmod(
+    cyan = mkinsub("SFORB", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO"), quiet = TRUE,
+    name = "sforb_path_1", dll_dir = "cyan_dlls", overwrite = TRUE),
+  hs_path_1 = mkinmod(
+    cyan = mkinsub("HS", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO"), quiet = TRUE,
+    name = "hs_path_1", dll_dir = "cyan_dlls", overwrite = TRUE)
+)
+cl_path_1 <- start_cluster(n_cores)

To obtain suitable starting values for the NLHM fits, separate pathway fits are performed for all datasets.

f_sep_1_const <- mmkin(
-  cyan_path_1,
-  cyan_ds,
-  error_model = "const",
-  cluster = cl_path_1,
-  quiet = TRUE)
-status(f_sep_1_const) |> kable()

	npar	AIC	BIC	Lik
SFO const

- - +

+f_sep_1_const <- mmkin(
+  cyan_path_1,
+  cyan_ds,
+  error_model = "const",
+  cluster = cl_path_1,
+  quiet = TRUE)
+status(f_sep_1_const) |> kable()

+ - - + @@ -1712,19 +1432,18 @@ status(f_sep_1_const) |> kable()

	Nambsheim	Tama	Gross-Umstadt	Sassafras	Lleida
sfo_path_1

f_sep_1_tc <- update(f_sep_1_const, error_model = "tc")
-status(f_sep_1_tc) |> kable()

- - - +

+f_sep_1_tc <- update(f_sep_1_const, error_model = "tc")
+status(f_sep_1_tc) |> kable()

+ - - + @@ -1778,18 +1497,18 @@ for the parent only fits is used as an argument no_random_effect to the mhmkin function. The possibility to do so was introduced in mkin version 1.2.2 which is currently under development.

f_saem_1 <- mhmkin(list(f_sep_1_const, f_sep_1_tc),
-  no_random_effect = illparms(cyan_saem_full),
-  cluster = cl_path_1)

status(f_saem_1) |> kable()

	Nambsheim	Tama	Gross-Umstadt	Sassafras	Lleida
sfo_path_1

- - +

+f_saem_1 <- mhmkin(list(f_sep_1_const, f_sep_1_tc),
+  no_random_effect = illparms(cyan_saem_full),
+  cluster = cl_path_1)

+status(f_saem_1) |> kable()

+ - - + @@ -1825,20 +1544,19 @@ Fisher Information Matrix could not be inverted for the fixed effects fits, ill-defined parameters cannot be determined using the illparms function, because it relies on the Fisher Information Matrix.

illparms(f_saem_1) |> kable()

	const	tc
sfo_path_1

+illparms(f_saem_1) |> kable()

---+++ - - + - - + @@ -1870,17 +1588,16 @@ sd(f_JCZ38_qlogis)

	const	tc
sfo_path_1

The model comparisons below suggest that the pathway fits using DFOP or SFORB for the parent compound provide the best fit.

anova(f_saem_1[, "const"]) |> kable(digits = 1)

- - - +

+anova(f_saem_1[, "const"]) |> kable(digits = 1)

+ - - + @@ -1919,17 +1636,16 @@ or SFORB for the parent compound provide the best fit.

	npar	AIC	BIC	Lik
sfo_path_1 const

anova(f_saem_1[1:4, ]) |> kable(digits = 1)

- - - +

+anova(f_saem_1[1:4, ]) |> kable(digits = 1)

+ - - + @@ -1991,27 +1707,29 @@ or SFORB for the parent compound provide the best fit.

	npar	AIC	BIC	Lik
sfo_path_1 const

For these two parent model, successful fits are shown below. Plots of the fits with the other parent models are shown in the Appendix.

plot(f_saem_1[["dfop_path_1", "tc"]])

+plot(f_saem_1[["dfop_path_1", "tc"]])

+ DFOP pathway fit with two-component error

DFOP pathway fit with two-component error

plot(f_saem_1[["sforb_path_1", "tc"]])

+plot(f_saem_1[["sforb_path_1", "tc"]])

+ SFORB pathway fit with two-component error

SFORB pathway fit with two-component error

A closer graphical analysis of these Figures shows that the residues of transformation product JCZ38 in the soils Tama and Nambsheim observed at later time points are strongly and systematically underestimated.

stopCluster(cl_path_1)

+stopCluster(cl_path_1)

Alternative pathway fits

Alternative pathway fits +

To improve the fit for JCZ38, a back-reaction from JSE76 to JCZ38 was introduced in an alternative version of the transformation pathway, in analogy to the back-reaction from K5A78 to K5A77. Both pairs of @@ -2020,56 +1738,55 @@ corresponding amide (Addendum 2014, p. 109). As FOMC provided the best fit for the parent, and the biexponential models DFOP and SFORB provided the best initial pathway fits, these three parent models are used in the alternative pathway fits.

cyan_path_2 <- list(
-  fomc_path_2 = mkinmod(
-    cyan = mkinsub("FOMC", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO", "JCZ38"),
-    name = "fomc_path_2", quiet = TRUE,
-    dll_dir = "cyan_dlls",
-    overwrite = TRUE
-  ),
-  dfop_path_2 = mkinmod(
-    cyan = mkinsub("DFOP", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO", "JCZ38"),
-    name = "dfop_path_2", quiet = TRUE,
-    dll_dir = "cyan_dlls",
-    overwrite = TRUE
-  ),
-  sforb_path_2 = mkinmod(
-    cyan = mkinsub("SFORB", c("JCZ38", "J9Z38")),
-    JCZ38 = mkinsub("SFO", "JSE76"),
-    J9Z38 = mkinsub("SFO"),
-    JSE76 = mkinsub("SFO", "JCZ38"),
-    name = "sforb_path_2", quiet = TRUE,
-    dll_dir = "cyan_dlls",
-    overwrite = TRUE
-  )
-)
-
-cl_path_2 <- start_cluster(n_cores)
-f_sep_2_const <- mmkin(
-  cyan_path_2,
-  cyan_ds,
-  error_model = "const",
-  cluster = cl_path_2,
-  quiet = TRUE)
-
-status(f_sep_2_const) |> kable()

- - - +

+cyan_path_2 <- list(
+  fomc_path_2 = mkinmod(
+    cyan = mkinsub("FOMC", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO", "JCZ38"),
+    name = "fomc_path_2", quiet = TRUE,
+    dll_dir = "cyan_dlls",
+    overwrite = TRUE
+  ),
+  dfop_path_2 = mkinmod(
+    cyan = mkinsub("DFOP", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO", "JCZ38"),
+    name = "dfop_path_2", quiet = TRUE,
+    dll_dir = "cyan_dlls",
+    overwrite = TRUE
+  ),
+  sforb_path_2 = mkinmod(
+    cyan = mkinsub("SFORB", c("JCZ38", "J9Z38")),
+    JCZ38 = mkinsub("SFO", "JSE76"),
+    J9Z38 = mkinsub("SFO"),
+    JSE76 = mkinsub("SFO", "JCZ38"),
+    name = "sforb_path_2", quiet = TRUE,
+    dll_dir = "cyan_dlls",
+    overwrite = TRUE
+  )
+)
+
+cl_path_2 <- start_cluster(n_cores)
+f_sep_2_const <- mmkin(
+  cyan_path_2,
+  cyan_ds,
+  error_model = "const",
+  cluster = cl_path_2,
+  quiet = TRUE)
+
+status(f_sep_2_const) |> kable()

+ - - + @@ -2099,19 +1816,18 @@ status(f_sep_2_const) |> kable()

	Nambsheim	Tama	Gross-Umstadt	Sassafras	Lleida
fomc_path_2

Using constant variance, separate fits converge with the exception of the fits to the Sassafras soil data.

f_sep_2_tc <- update(f_sep_2_const, error_model = "tc")
-status(f_sep_2_tc) |> kable()

- - - +

+f_sep_2_tc <- update(f_sep_2_const, error_model = "tc")
+status(f_sep_2_tc) |> kable()

+ - - + @@ -2142,18 +1858,18 @@ status(f_sep_2_tc) |> kable()

Using the two-component error model, all separate fits converge with the exception of the alternative pathway fit with DFOP used for the parent and the Sassafras dataset.

f_saem_2 <- mhmkin(list(f_sep_2_const, f_sep_2_tc),
-  no_random_effect = illparms(cyan_saem_full[2:4, ]),
-  cluster = cl_path_2)

status(f_saem_2) |> kable()

	Nambsheim	Tama	Gross-Umstadt	Sassafras	Lleida
fomc_path_2

- - +

+f_saem_2 <- mhmkin(list(f_sep_2_const, f_sep_2_tc),
+  no_random_effect = illparms(cyan_saem_full[2:4, ]),
+  cluster = cl_path_2)

+status(f_saem_2) |> kable()

+ - - + @@ -2175,20 +1891,19 @@ parent and the Sassafras dataset.

The hierarchical fits for the alternative pathway completed successfully, with the exception of the model using FOMC for the parent compound and constant variance as the error model.

illparms(f_saem_2) |> kable()

	const	tc
fomc_path_2

+illparms(f_saem_2) |> kable()

---+++ - - + - - + @@ -2211,17 +1926,16 @@ compound and constant variance as the error model.

random effects for the formation fractions for the pathways from JCZ38 to JSE76, and for the reverse pathway from JSE76 to JCZ38 are ill-defined.

anova(f_saem_2[, "tc"]) |> kable(digits = 1)

	const	tc
fomc_path_2

- - +

+anova(f_saem_2[, "tc"]) |> kable(digits = 1)

+ - - + @@ -2246,17 +1960,16 @@ ill-defined.

	npar	AIC	BIC	Lik
fomc_path_2 tc

anova(f_saem_2[2:3,]) |> kable(digits = 1)

- - - +

+anova(f_saem_2[2:3,]) |> kable(digits = 1)

+ - - + @@ -2294,30 +2007,31 @@ and BIC values and are plotted below. Compared with the original pathway, the AIC and BIC values indicate a large improvement. This is confirmed by the plots, which show that the metabolite JCZ38 is fitted much better with this model.

plot(f_saem_2[["fomc_path_2", "tc"]])

+plot(f_saem_2[["fomc_path_2", "tc"]])

+ FOMC pathway fit with two-component error, alternative pathway

FOMC pathway fit with two-component error, alternative pathway

plot(f_saem_2[["dfop_path_2", "tc"]])

+plot(f_saem_2[["dfop_path_2", "tc"]])

+ DFOP pathway fit with two-component error, alternative pathway

DFOP pathway fit with two-component error, alternative pathway

plot(f_saem_2[["sforb_path_2", "tc"]])

+plot(f_saem_2[["sforb_path_2", "tc"]])

+ SFORB pathway fit with two-component error, alternative pathway

SFORB pathway fit with two-component error, alternative pathway

Refinement of alternative pathway fits

All ill-defined random effects that were identified in the parent only fits and in the above pathway fits, are excluded for the final evaluations below. For this purpose, a list of character vectors is @@ -2325,29 +2039,29 @@ created below that can be indexed by row and column indices, and which contains the degradation parameter names for which random effects should be excluded for each of the hierarchical fits contained in f_saem_2.

no_ranef <- matrix(list(), nrow = 3, ncol = 2, dimnames = dimnames(f_saem_2))
-no_ranef[["fomc_path_2", "const"]] <- c("log_beta", "f_JCZ38_qlogis", "f_JSE76_qlogis")
-no_ranef[["fomc_path_2", "tc"]] <- c("cyan_0", "f_JCZ38_qlogis", "f_JSE76_qlogis")
-no_ranef[["dfop_path_2", "const"]] <- c("cyan_0", "f_JCZ38_qlogis", "f_JSE76_qlogis")
-no_ranef[["dfop_path_2", "tc"]] <- c("cyan_0", "log_k1", "f_JCZ38_qlogis", "f_JSE76_qlogis")
-no_ranef[["sforb_path_2", "const"]] <- c("cyan_free_0",
-  "f_JCZ38_qlogis", "f_JSE76_qlogis")
-no_ranef[["sforb_path_2", "tc"]] <- c("cyan_free_0", "log_k_cyan_free_bound",
-  "f_JCZ38_qlogis", "f_JSE76_qlogis")
-clusterExport(cl_path_2, "no_ranef")
-
-f_saem_3 <- update(f_saem_2,
-  no_random_effect = no_ranef,
-  cluster = cl_path_2)

status(f_saem_3) |> kable()

	npar	AIC	BIC	Lik
dfop_path_2 const

- - +

+no_ranef <- matrix(list(), nrow = 3, ncol = 2, dimnames = dimnames(f_saem_2))
+no_ranef[["fomc_path_2", "const"]] <- c("log_beta", "f_JCZ38_qlogis", "f_JSE76_qlogis")
+no_ranef[["fomc_path_2", "tc"]] <- c("cyan_0", "f_JCZ38_qlogis", "f_JSE76_qlogis")
+no_ranef[["dfop_path_2", "const"]] <- c("cyan_0", "f_JCZ38_qlogis", "f_JSE76_qlogis")
+no_ranef[["dfop_path_2", "tc"]] <- c("cyan_0", "log_k1", "f_JCZ38_qlogis", "f_JSE76_qlogis")
+no_ranef[["sforb_path_2", "const"]] <- c("cyan_free_0",
+  "f_JCZ38_qlogis", "f_JSE76_qlogis")
+no_ranef[["sforb_path_2", "tc"]] <- c("cyan_free_0", "log_k_cyan_free_bound",
+  "f_JCZ38_qlogis", "f_JSE76_qlogis")
+clusterExport(cl_path_2, "no_ranef")
+
+f_saem_3 <- update(f_saem_2,
+  no_random_effect = no_ranef,
+  cluster = cl_path_2)

+status(f_saem_3) |> kable()

+ - - + @@ -2370,15 +2084,14 @@ f_saem_3 <- update(f_saem_2, all updated fits completed successfully. However, the Fisher Information Matrix for the fixed effects (Fth) could not be inverted, so no confidence intervals for the optimised parameters are available.

illparms(f_saem_3) |> kable()

	const	tc
fomc_path_2

- - +

+illparms(f_saem_3) |> kable()

+ - - + @@ -2397,17 +2110,16 @@ confidence intervals for the optimised parameters are available.

	const	tc
fomc_path_2

anova(f_saem_3[, "tc"]) |> kable(digits = 1)

- - - +

+anova(f_saem_3[, "tc"]) |> kable(digits = 1)

+ - - + @@ -2432,17 +2144,16 @@ confidence intervals for the optimised parameters are available.

	npar	AIC	BIC	Lik
fomc_path_2 tc

anova(f_saem_3[2:3,]) |> kable(digits = 1)

- - - +

+anova(f_saem_3[2:3,]) |> kable(digits = 1)

+ - - + @@ -2478,11 +2189,13 @@ confidence intervals for the optimised parameters are available.

two-component error) are lower than in the previous fits with the alternative pathway, the practical value of these refined evaluations is limited as no confidence intervals are obtained.

stopCluster(cl_path_2)

+stopCluster(cl_path_2)

Conclusion

Conclusion +

It was demonstrated that a relatively complex transformation pathway with parallel formation of two primary metabolites and one secondary metabolite can be fitted even if the data in the individual datasets are @@ -2492,50 +2205,55 @@ practical feasibility of iterative refinements based on ill-defined parameters and of alternative checks of parameter identifiability based on multistart runs.

Acknowledgements

Acknowledgements +

The helpful comments by Janina Wöltjen of the German Environment Agency are gratefully acknowledged.

Appendix

Plots of fits that were not refined further

plot(f_saem_1[["sfo_path_1", "tc"]])

Appendix +

Plots of fits that were not refined further +

+plot(f_saem_1[["sfo_path_1", "tc"]])

+ SFO pathway fit with two-component error

SFO pathway fit with two-component error

plot(f_saem_1[["fomc_path_1", "tc"]])

+plot(f_saem_1[["fomc_path_1", "tc"]])

+ FOMC pathway fit with two-component error

FOMC pathway fit with two-component error

plot(f_saem_1[["sforb_path_1", "tc"]])

+plot(f_saem_1[["sforb_path_1", "tc"]])

+ HS pathway fit with two-component error

HS pathway fit with two-component error

Hierarchical fit listings

Pathway 1

Hierarchical fit listings +

Pathway 1 +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:31:33 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:35:03 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - k_cyan * cyan
@@ -2548,7 +2266,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 480.873 s
+Fitted in 445.08 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -2653,17 +2371,17 @@ JCZ38  14.46  48.05
 J9Z38 103.86 345.00
 JSE76 143.91 478.04
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:32:28 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:34:47 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - k_cyan * cyan
@@ -2676,7 +2394,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 534.75 s
+Fitted in 428.46 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -2783,17 +2501,17 @@ JCZ38  15.81  52.52
 J9Z38 107.97 358.68
 JSE76 114.20 379.35
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:33:51 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:35:21 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - (alpha/beta) * 1/((time/beta) + 1) * cyan
@@ -2808,7 +2526,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 618.676 s
+Fitted in 462.739 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -2928,17 +2646,17 @@ JCZ38  20.53   68.19       NA
 J9Z38 140.07  465.32       NA
 JSE76 318.86 1059.22       NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:34:01 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:35:31 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - (alpha/beta) * 1/((time/beta) + 1) * cyan
@@ -2953,7 +2671,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 627.822 s
+Fitted in 472.067 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -3066,17 +2784,17 @@ JCZ38  22.01  73.1       NA
 J9Z38 130.09 432.2       NA
 JSE76 210.98 700.9       NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:33:18 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:36:02 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -3095,7 +2813,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 584.724 s
+Fitted in 503.671 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -3212,17 +2930,17 @@ JCZ38  22.34  74.22       NA      NA      NA
 J9Z38 119.92 398.36       NA      NA      NA
 JSE76 200.41 665.76       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:35:43 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:38:54 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -3241,7 +2959,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 729.575 s
+Fitted in 674.81 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -3358,17 +3076,17 @@ JCZ38  22.20  73.75       NA      NA      NA
 J9Z38 108.23 359.53       NA      NA      NA
 JSE76 179.30 595.62       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:34:05 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:36:26 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -3386,7 +3104,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 632.71 s
+Fitted in 527.639 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -3523,17 +3241,17 @@ JCZ38  22.98  76.33       NA           NA           NA
 J9Z38 116.28 386.29       NA           NA           NA
 JSE76 193.42 642.53       NA           NA           NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:37:01 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:38:31 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -3551,7 +3269,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 807.852 s
+Fitted in 651.67 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -3688,17 +3406,17 @@ JCZ38  22.34  74.2       NA           NA           NA
 J9Z38 110.14 365.9       NA           NA           NA
 JSE76 177.11 588.3       NA           NA           NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:33:29 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:36:13 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ifelse(time <= tb, k1, k2) * cyan
@@ -3713,7 +3431,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 596.235 s
+Fitted in 514.353 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -3830,20 +3548,21 @@ JCZ38  22.27  73.98       NA      NA      NA
 J9Z38 113.09 375.69       NA      NA      NA
 JSE76 187.01 621.23       NA      NA      NA
 
-

+ +

Pathway 2

Pathway 2 +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:46:08 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:47:30 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - (alpha/beta) * 1/((time/beta) + 1) * cyan
@@ -3858,7 +3577,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 536.687 s
+Fitted in 505.533 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -3999,17 +3718,17 @@ JCZ38   7.075  23.50       NA
 J9Z38 117.249 389.49       NA
 JSE76  14.169  47.07       NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:47:06 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:48:20 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -4029,7 +3748,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 594.209 s
+Fitted in 555.413 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -4190,17 +3909,17 @@ JCZ38  12.50  41.53       NA      NA      NA
 J9Z38 118.69 394.27       NA      NA      NA
 JSE76  24.32  80.78       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:49:43 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:51:02 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -4220,7 +3939,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 751.883 s
+Fitted in 717.675 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -4381,17 +4100,17 @@ JCZ38   8.93  29.66       NA      NA      NA
 J9Z38 110.45 366.89       NA      NA      NA
 JSE76  17.96  59.66       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:46:57 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:48:15 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -4409,7 +4128,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 585.771 s
+Fitted in 550.772 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -4577,17 +4296,17 @@ JCZ38   7.203  23.93       NA           NA           NA
 J9Z38 131.918 438.22       NA           NA           NA
 JSE76  14.287  47.46       NA           NA           NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 18:50:00 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 20:50:59 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -4605,7 +4324,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 767.874 s
+Fitted in 714.48 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -4773,21 +4492,22 @@ JCZ38   8.535  28.35       NA           NA           NA
 J9Z38 105.517 350.52       NA           NA           NA
 JSE76  17.837  59.25       NA           NA           NA
 
-

+ +

Pathway 2, refined fits

Pathway 2, refined fits +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 19:03:52 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:04:10 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - (alpha/beta) * 1/((time/beta) + 1) * cyan
@@ -4802,7 +4522,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 830.375 s
+Fitted in 786.038 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -4919,18 +4639,18 @@ JCZ38  12.03  39.96       NA
 J9Z38 111.14 369.19       NA
 JSE76  23.77  78.98       NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 19:05:47 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:05:56 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -4950,7 +4670,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 945.728 s
+Fitted in 892.139 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -5085,18 +4805,18 @@ JCZ38  13.04  43.33       NA      NA      NA
 J9Z38 120.93 401.73       NA      NA      NA
 JSE76  26.39  87.68       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 19:05:49 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:02 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -5116,7 +4836,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 947.743 s
+Fitted in 898.534 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -5251,18 +4971,18 @@ JCZ38  11.49  38.18       NA      NA      NA
 J9Z38 107.55 357.28       NA      NA      NA
 JSE76  27.20  90.36       NA      NA      NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 19:05:38 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:05:52 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -5280,7 +5000,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 936.368 s
+Fitted in 888.333 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Constant variance 
@@ -5422,18 +5142,18 @@ JCZ38  12.92  42.93       NA           NA           NA
 J9Z38 114.71 381.07       NA           NA           NA
 JSE76  26.04  86.51       NA           NA           NA
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 19:05:52 2025 
-Date of summary: Thu Feb 13 19:05:53 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:21 2025 
+Date of summary: Mon May 12 21:06:22 2025 
 
 Equations:
 d_cyan_free/dt = - k_cyan_free * cyan_free - k_cyan_free_bound *
@@ -5451,7 +5171,7 @@ Data:
 
 Model predictions using solution type deSolve 
 
-Fitted in 950.661 s
+Fitted in 916.619 s
 Using 300, 100 iterations and 10 chains
 
 Variance model: Two-component variance function 
@@ -5593,19 +5313,20 @@ JCZ38  11.06  36.75       NA           NA           NA
 J9Z38 106.71 354.49       NA           NA           NA
 JSE76  25.44  84.51       NA           NA           NA
 
-

+ +

Session info

R version 4.4.2 (2024-10-31)
+
+Session info
+
+R version 4.5.0 (2025-04-11)
 Platform: x86_64-pc-linux-gnu
 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
+LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0  LAPACK version 3.11.0
 
 locale:
  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
@@ -5623,75 +5344,52 @@ attached base packages:
 [8] base     
 
 other attached packages:
-[1] saemix_3.3      npde_3.5        knitr_1.49      mkin_1.2.9     
-[5] rmarkdown_2.29  nvimcom_0.9-167
+[1] rmarkdown_2.29  nvimcom_0.9-167 saemix_3.3      npde_3.5       
+[5] knitr_1.49      mkin_1.2.10    
 
 loaded via a namespace (and not attached):
- [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    lattice_0.22-6   
- [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.1    grid_4.4.2       
- [9] fastmap_1.2.0     cellranger_1.1.0  jsonlite_1.8.9    processx_3.8.4   
-[13] pkgbuild_1.4.5    deSolve_1.40      mclust_6.1.1      ps_1.8.1         
-[17] gridExtra_2.3     fansi_1.0.6       scales_1.3.0      codetools_0.2-20 
-[21] jquerylib_0.1.4   cli_3.6.3         rlang_1.1.4       munsell_0.5.1    
-[25] cachem_1.1.0      yaml_2.3.10       inline_0.3.20     tools_4.4.2      
-[29] colorout_1.3-2    dplyr_1.1.4       colorspace_2.1-1  ggplot2_3.5.1    
-[33] vctrs_0.6.5       R6_2.5.1          zoo_1.8-12        lifecycle_1.0.4  
-[37] MASS_7.3-61       pkgconfig_2.0.3   callr_3.7.6       pillar_1.9.0     
-[41] bslib_0.8.0       gtable_0.3.6      glue_1.8.0        xfun_0.49        
-[45] tibble_3.2.1      lmtest_0.9-40     tidyselect_1.2.1  htmltools_0.5.8.1
-[49] nlme_3.1-166      compiler_4.4.2    readxl_1.4.3     
+ [1] sass_0.4.9        generics_0.1.3    lattice_0.22-6    digest_0.6.37    
+ [5] magrittr_2.0.3    evaluate_1.0.3    grid_4.5.0        fastmap_1.2.0    
+ [9] cellranger_1.1.0  jsonlite_1.9.0    processx_3.8.6    pkgbuild_1.4.6   
+[13] deSolve_1.40      mclust_6.1.1      ps_1.9.0          gridExtra_2.3    
+[17] scales_1.3.0      codetools_0.2-20  textshaping_1.0.0 jquerylib_0.1.4  
+[21] cli_3.6.4         rlang_1.1.5       munsell_0.5.1     cachem_1.1.0     
+[25] yaml_2.3.10       inline_0.3.21     tools_4.5.0       dplyr_1.1.4      
+[29] colorspace_2.1-1  ggplot2_3.5.1     vctrs_0.6.5       R6_2.6.1         
+[33] zoo_1.8-13        lifecycle_1.0.4   fs_1.6.5          htmlwidgets_1.6.4
+[37] MASS_7.3-65       ragg_1.3.3        pkgconfig_2.0.3   desc_1.4.3       
+[41] callr_3.7.6       pkgdown_2.1.1     pillar_1.10.1     bslib_0.9.0      
+[45] gtable_0.3.6      glue_1.8.0        systemfonts_1.2.1 xfun_0.51        
+[49] tibble_3.2.1      lmtest_0.9-40     tidyselect_1.2.1  htmltools_0.5.8.1
+[53] nlme_3.1-168      compiler_4.5.0    readxl_1.4.4

Hardware info

Hardware info +

CPU model: AMD Ryzen 9 7950X 16-Core Processor

MemTotal:       64927788 kB

MemTotal:       64927780 kB

- - - - + - +

- + - - - - - + diff --git a/docs/articles/prebuilt/2022_dmta_parent.html b/docs/articles/prebuilt/2022_dmta_parent.html index 14f12c0e..10035c99 100644 --- a/docs/articles/prebuilt/2022_dmta_parent.html +++ b/docs/articles/prebuilt/2022_dmta_parent.html @@ -1,396 +1,100 @@ - - - + - - - - - - - - - -Testing hierarchical parent degradation kinetics with residue data on dimethenamid and dimethenamid-P - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Testing hierarchical parent degradation kinetics with residue data on dimethenamid and dimethenamid-P • mkin + + + + + - + + + +

Testing hierarchical parent degradation -kinetics with residue data on dimethenamid and dimethenamid-P

Johannes Ranke

Last change on 5 January 2023, last compiled on 13 -Februar 2025

Testing hierarchical parent degradation kinetics with residue data on dimethenamid and dimethenamid-P

Johannes +Ranke

+ +

Last change on 5 January +2023, last compiled on 12 Mai 2025

+ + Source: vignettes/prebuilt/2022_dmta_parent.rmd +

2022_dmta_parent.rmd

- -

-true -

- -

Introduction

+ + +

Introduction +

The purpose of this document is to demonstrate how nonlinear hierarchical models (NLHM) based on the parent degradation models SFO, FOMC, DFOP and HS can be fitted with the mkin package.

@@ -398,7 +102,7 @@ FOMC, DFOP and HS can be fitted with the mkin package.

173340 (Application of nonlinear hierarchical models to the kinetic evaluation of chemical degradation data) of the German Environment Agency carried out in 2022 and 2023.

The mkin package is used in version 1.2.9. It contains the test data +

The mkin package is used in version 1.2.10. It contains the test data and the functions used in 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.

@@ -406,19 +110,21 @@ make the convergence plot function available.

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)
-}

+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)
+}

Data

Data +

The test data are available in the mkin package as an object of class mkindsg (mkin dataset group) under the identifier dimethenamid_2018. The following preprocessing steps are @@ -440,37 +146,37 @@ and Elliot 2) are combined, resulting in dimethenamid (DMTA) data from six soils.

The following commented R code performs this preprocessing.

# Apply a function to each of the seven datasets in the mkindsg object to create a list
-dmta_ds <- lapply(1:7, function(i) {
-  ds_i <- dimethenamid_2018$ds[[i]]$data                     # Get a dataset
-  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"              # Rename DMTAP to DMTA
-  ds_i <- subset(ds_i, name == "DMTA", c("name", "time", "value")) # Select data
-  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]  # Normalise time
-  ds_i                                                       # Return the dataset
-})
-
-# Use dataset titles as names for the list elements
-names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-
-# Combine data for Elliot soil to obtain a named list with six elements
-dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) #
-dmta_ds[["Elliot 1"]] <- NULL
-dmta_ds[["Elliot 2"]] <- NULL

+# Apply a function to each of the seven datasets in the mkindsg object to create a list
+dmta_ds <- lapply(1:7, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data                     # Get a dataset
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"              # Rename DMTAP to DMTA
+  ds_i <- subset(ds_i, name == "DMTA", c("name", "time", "value")) # Select data
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]  # Normalise time
+  ds_i                                                       # Return the dataset
+})
+
+# Use dataset titles as names for the list elements
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+
+# Combine data for Elliot soil to obtain a named list with six elements
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) #
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL

The following tables show the 6 datasets.

for (ds_name in names(dmta_ds)) {
-    print(kable(mkin_long_to_wide(dmta_ds[[ds_name]]),
-      caption = paste("Dataset", ds_name),
-      label = paste0("tab:", ds_name), booktabs = TRUE))
-    cat("\n\\clearpage\n")
-}

Hierarchical SFO path 1 fit with constant variance
	npar	AIC	BIC	Lik
dfop_path_2 const

+for (ds_name in names(dmta_ds)) {
+    print(kable(mkin_long_to_wide(dmta_ds[[ds_name]]),
+      caption = paste("Dataset", ds_name),
+      label = paste0("tab:", ds_name), booktabs = TRUE))
+    cat("\n\\clearpage\n")
+}

- - + - - + @@ -502,14 +208,12 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Calke
time	DMTA
0

- +

- - + - - + @@ -577,14 +281,12 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Borstel
time	DMTA
0.000000

- +

- - + - - + @@ -732,14 +434,12 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Flaach
time	DMTA
0.0000000

- +

- - + - - + @@ -839,14 +539,12 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset BBA 2.2
time	DMTA
0.0000000

- +

- - + - - + @@ -946,14 +644,12 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset BBA 2.3
time	DMTA
0.0000000

- +

- - + - - + @@ -1150,24 +846,25 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Elliot
time	DMTA
0.000000

Separate evaluations

Separate evaluations +

In order to obtain suitable starting parameters for the NLHM fits, separate fits of the four 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", "HS")
-f_sep_const <- mmkin(
-  deg_mods,
-  dmta_ds,
-  error_model = "const",
-  quiet = TRUE)
-
-status(f_sep_const) |> kable()

- - - +

+deg_mods <- c("SFO", "FOMC", "DFOP", "HS")
+f_sep_const <- mmkin(
+  deg_mods,
+  dmta_ds,
+  error_model = "const",
+  quiet = TRUE)
+
+status(f_sep_const) |> kable()

+ @@ -1175,8 +872,7 @@ status(f_sep_const) |> kable() - - + @@ -1221,11 +917,11 @@ to converge. All separate fits with constant variance converged, with the sole exception of the HS fit to the BBA 2.2 data. To prepare for fitting NLHM using the two-component error model, the separate fits are updated assuming two-component error.

f_sep_tc <- update(f_sep_const, error_model = "tc")
-status(f_sep_tc) |> kable()

	Calke	Borstel	BBA 2.2	BBA 2.3	Elliot
SFO

- - +

+f_sep_tc <- update(f_sep_const, error_model = "tc")
+status(f_sep_tc) |> kable()

+ @@ -1233,8 +929,7 @@ status(f_sep_tc) |> kable() - - + @@ -1278,8 +973,9 @@ status(f_sep_tc) |> kable() converge with constant variance did converge, but other non-SFO fits failed to converge.

Hierarchichal model fits

Hierarchichal model fits +

The following code fits eight versions of hierarchical models to the data, using SFO, FOMC, DFOP and HS for the parent compound, and using either constant variance or two-component error for the error model. The @@ -1291,18 +987,18 @@ all eight versions in parallel (given a sufficient number of computing cores being available) to save execution time.

Convergence plots and summaries for these fits are shown in the appendix.

f_saem <- mhmkin(list(f_sep_const, f_sep_tc), transformations = "saemix")

+f_saem <- mhmkin(list(f_sep_const, f_sep_tc), transformations = "saemix")

The output of the status function shows that all fits terminated successfully.

status(f_saem) |> kable()

	Calke	Borstel	BBA 2.2	BBA 2.3	Elliot
SFO

- - +

+status(f_saem) |> kable()

+ - - + @@ -1328,17 +1024,16 @@ terminated successfully.

	const	tc
SFO

The AIC and BIC values show that the biphasic models DFOP and HS give the best fits.

anova(f_saem) |> kable(digits = 1)

- - - +

+anova(f_saem) |> kable(digits = 1)

+ - - + @@ -1404,22 +1099,22 @@ it shows the lowest AIC and BIC values in the first set of fits when combined with the two-component error model. Therefore, the DFOP model was selected for further refinements of the fits with the aim to make the model fully identifiable.

Parameter identifiability based on the Fisher Information -Matrix

Parameter identifiability based on the Fisher Information +Matrix +

Using the illparms function, ill-defined statistical model parameters such as standard deviations of the degradation parameters in the population and error model parameters can be found.

illparms(f_saem) |> kable()

	npar	AIC	BIC	Lik
SFO const

- - +

+illparms(f_saem) |> kable()

+ - - + @@ -1450,32 +1145,34 @@ be obtained for this standard deviation when using different starting values.

The thus identified overparameterisation is addressed by removing the random effect for k2 from the parameter model.

f_saem_dfop_tc_no_ranef_k2 <- update(f_saem[["DFOP", "tc"]],
-  no_random_effect = "k2")

+f_saem_dfop_tc_no_ranef_k2 <- update(f_saem[["DFOP", "tc"]],
+  no_random_effect = "k2")

For the resulting fit, it is checked whether there are still ill-defined parameters,

illparms(f_saem_dfop_tc_no_ranef_k2)

+illparms(f_saem_dfop_tc_no_ranef_k2)

which is not the case. Below, the refined model is compared with the previous best model. The model without random effect for k2 is a reduced version of the previous model. Therefore, the models are nested and can be compared using the likelihood ratio test. This is achieved with the argument test = TRUE to the anova function.

anova(f_saem[["DFOP", "tc"]], f_saem_dfop_tc_no_ranef_k2, test = TRUE) |>
-  kable(format.args = list(digits = 4))

	const	tc
SFO

+anova(f_saem[["DFOP", "tc"]], f_saem_dfop_tc_no_ranef_k2, test = TRUE) |>
+  kable(format.args = list(digits = 4))

--------++++++++ - - + @@ -1484,8 +1181,7 @@ achieved with the argument test = TRUE to the - - + @@ -1513,35 +1209,35 @@ achieved with the argument test = TRUE to the random effect for k2. The p value of the likelihood ratio test is much greater than 0.05, indicating that the model with the higher likelihood (here the model with random effects for all -degradation parameters f_saem[["DFOP", "tc"]]) does not fit +degradation parameters f_saem[["DFOP", "tc"]]) does not fit significantly better than the model with the lower likelihood (the reduced model f_saem_dfop_tc_no_ranef_k2).

Therefore, AIC, BIC and likelihood ratio test suggest the use of the reduced model.

The convergence of the fit is checked visually.

+ Convergence plot for the NLHM DFOP fit with two-component error and without a random effect on 'k2'

Convergence plot for the NLHM DFOP fit with two-component error and without a random effect on ‘k2’

All parameters appear to have converged to a satisfactory degree. The final fit is plotted using the plot method from the mkin package.

plot(f_saem_dfop_tc_no_ranef_k2)

+plot(f_saem_dfop_tc_no_ranef_k2)

+ Plot of the final NLHM DFOP fit

Plot of the final NLHM DFOP fit

Finally, a summary report of the fit is produced.

summary(f_saem_dfop_tc_no_ranef_k2)

+summary(f_saem_dfop_tc_no_ranef_k2)

saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:33 2025 
-Date of summary: Thu Feb 13 16:33:34 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:46 2025 
+Date of summary: Mon May 12 21:06:46 2025 
 
 Equations:
 d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1553,7 +1249,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 3.778 s
+Fitted in 3.953 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
@@ -1615,8 +1311,9 @@ Estimated disappearance times:
       DT50  DT90 DT50back DT50_k1 DT50_k2
 DMTA 11.45 41.32    12.44   10.82   81.85

Alternative check of parameter identifiability

Alternative check of parameter identifiability +

The parameter check used in the illparms function is based on a quadratic approximation of the likelihood surface near its optimum, which is calculated using the Fisher Information Matrix (FIM). @@ -1626,12 +1323,13 @@ approach has recently been implemented in mkin.

of the saem algorithm with different parameter combinations, sampled from the range of the parameters obtained for the individual datasets fitted separately using nonlinear regression.

f_saem_dfop_tc_multi <- multistart(f_saem[["DFOP", "tc"]], n = 50, cores = 15)

par(mar = c(6.1, 4.1, 2.1, 2.1))
-parplot(f_saem_dfop_tc_multi, lpos = "bottomright", ylim = c(0.3, 10), las = 2)

+f_saem_dfop_tc_multi <- multistart(f_saem[["DFOP", "tc"]], n = 50, cores = 15)

+par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_multi, lpos = "bottomright", ylim = c(0.3, 10), las = 2)

+ Scaled parameters from the multistart runs, full model

Scaled parameters from the multistart runs, full model

@@ -1642,34 +1340,36 @@ parameter g.

The parameter boxplots of the multistart runs with the reduced model shown below indicate that all runs give similar results, regardless of the starting parameters.

f_saem_dfop_tc_no_ranef_k2_multi <- multistart(f_saem_dfop_tc_no_ranef_k2,
-  n = 50, cores = 15)

par(mar = c(6.1, 4.1, 2.1, 2.1))
-parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2,
-  lpos = "bottomright")

+f_saem_dfop_tc_no_ranef_k2_multi <- multistart(f_saem_dfop_tc_no_ranef_k2,
+  n = 50, cores = 15)

+par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2,
+  lpos = "bottomright")

+ Scaled parameters from the multistart runs, reduced model

Scaled parameters from the multistart runs, reduced model

When only the parameters of the top 25% of the fits are shown (based on a feature introduced in mkin 1.2.2 currently under development), the scatter is even less as shown below.

par(mar = c(6.1, 4.1, 2.1, 2.1))
-parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2, llquant = 0.25,
-  lpos = "bottomright")

+par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2, llquant = 0.25,
+  lpos = "bottomright")

- $Scaled parameters from the multistart runs, reduced model, fits with the top 25\% likelihood values$ -

+ $Scaled parameters from the multistart runs, reduced model, fits with the top 25\% likelihood values$

Scaled parameters from the multistart runs, reduced model, fits with the top 25% likelihood values

Conclusions

Conclusions +

Fitting the four parent degradation models SFO, FOMC, DFOP and HS as part of hierarchical model fits with two different error models and normal distributions of the transformed degradation parameters works @@ -1681,35 +1381,39 @@ was fully identifiable and showed the lowest values for the model selection criteria AIC and BIC. The reliability of the identification of all model parameters was confirmed using multiple starting values.

Acknowledgements

Acknowledgements +

The helpful comments by Janina Wöltjen of the German Environment Agency are gratefully acknowledged.

References

References +

Duchesne, Ronan, Anissa Guillemin, Olivier Gandrillon, and Fabien Crauste. 2021. “Practical Identifiability in the Frame of Nonlinear Mixed Effects Models: The Example of the in Vitro -Erythropoiesis.” BMC Bioinformatics 22 (478). https://doi.org/10.1186/s12859-021-04373-4. +Erythropoiesis.” BMC Bioinformatics 22 (478). https://doi.org/10.1186/s12859-021-04373-4.

Appendix

Hierarchical model fit listings

Appendix +

Hierarchical model fit listings +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:26 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:38 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - k_DMTA * DMTA
@@ -1719,7 +1423,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 0.792 s
+Fitted in 0.793 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Constant variance 
@@ -1771,17 +1475,17 @@ Estimated disappearance times:
       DT50  DT90
 DMTA 12.24 40.65
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:28 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:40 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - k_DMTA * DMTA
@@ -1791,7 +1495,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 2.245 s
+Fitted in 2.462 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
@@ -1845,17 +1549,17 @@ Estimated disappearance times:
       DT50  DT90
 DMTA 12.24 40.65
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:27 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:39 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
@@ -1865,7 +1569,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 1.409 s
+Fitted in 1.408 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Constant variance 
@@ -1922,17 +1626,17 @@ Estimated disappearance times:
       DT50  DT90 DT50back
 DMTA 11.41 42.53     12.8
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:28 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:40 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
@@ -1942,7 +1646,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 2.811 s
+Fitted in 2.735 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
@@ -2001,17 +1705,17 @@ Estimated disappearance times:
       DT50  DT90 DT50back
 DMTA 11.05 42.81    12.89
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:27 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:39 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -2023,7 +1727,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 1.638 s
+Fitted in 1.769 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Constant variance 
@@ -2085,17 +1789,17 @@ Estimated disappearance times:
       DT50 DT90 DT50back DT50_k1 DT50_k2
 DMTA 11.79 42.8    12.88   11.09   76.46
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:28 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:41 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -2107,7 +1811,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 3.024 s
+Fitted in 3.069 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
@@ -2171,17 +1875,17 @@ Estimated disappearance times:
       DT50  DT90 DT50back DT50_k1 DT50_k2
 DMTA 11.39 41.36    12.45   10.74   83.48
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:28 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:39 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
@@ -2191,7 +1895,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 2.301 s
+Fitted in 1.996 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Constant variance 
@@ -2253,17 +1957,17 @@ Estimated disappearance times:
       DT50  DT90 DT50back DT50_k1 DT50_k2
 DMTA 8.801 67.91    20.44   8.801   31.13
 
-

+ +


 saemix version used for fitting:      3.3 
-mkin version used for pre-fitting:  1.2.9 
-R version used for fitting:         4.4.2 
-Date of fit:     Thu Feb 13 16:33:29 2025 
-Date of summary: Thu Feb 13 16:34:39 2025 
+mkin version used for pre-fitting:  1.2.10 
+R version used for fitting:         4.5.0 
+Date of fit:     Mon May 12 21:06:41 2025 
+Date of summary: Mon May 12 21:07:51 2025 
 
 Equations:
 d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
@@ -2273,7 +1977,7 @@ Data:
 
 Model predictions using solution type analytical 
 
-Fitted in 3.264 s
+Fitted in 3.389 s
 Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
@@ -2337,69 +2041,63 @@ Estimated disappearance times:
       DT50  DT90 DT50back DT50_k1 DT50_k2
 DMTA 11.84 51.71    15.57   11.84   29.66
 
-

+ +

Hierarchical model convergence plots

Hierarchical model convergence plots +

+ Convergence plot for the NLHM SFO fit with constant variance

Convergence plot for the NLHM SFO fit with constant variance

+ Convergence plot for the NLHM SFO fit with two-component error

Convergence plot for the NLHM SFO fit with two-component error

+ Convergence plot for the NLHM FOMC fit with constant variance

Convergence plot for the NLHM FOMC fit with constant variance

+ Convergence plot for the NLHM FOMC fit with two-component error

Convergence plot for the NLHM FOMC fit with two-component error

+ Convergence plot for the NLHM DFOP fit with constant variance

Convergence plot for the NLHM DFOP fit with constant variance

+ Convergence plot for the NLHM DFOP fit with two-component error

Convergence plot for the NLHM DFOP fit with two-component error

+ Convergence plot for the NLHM HS fit with constant variance

Convergence plot for the NLHM HS fit with constant variance

+ Convergence plot for the NLHM HS fit with two-component error

Convergence plot for the NLHM HS fit with two-component error

Session info

R version 4.4.2 (2024-10-31)
+
+Session info
+
+R version 4.5.0 (2025-04-11)
 Platform: x86_64-pc-linux-gnu
 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
+LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0  LAPACK version 3.11.0
 
 locale:
  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
@@ -2417,73 +2115,50 @@ attached base packages:
 [8] base     
 
 other attached packages:
-[1] saemix_3.3      npde_3.5        knitr_1.49      mkin_1.2.9     
-[5] rmarkdown_2.29  nvimcom_0.9-167
+[1] rmarkdown_2.29  nvimcom_0.9-167 saemix_3.3      npde_3.5       
+[5] knitr_1.49      mkin_1.2.10    
 
 loaded via a namespace (and not attached):
- [1] jsonlite_1.8.9    gtable_0.3.6      dplyr_1.1.4       compiler_4.4.2   
- [5] tidyselect_1.2.1  colorout_1.3-2    tinytex_0.54      gridExtra_2.3    
- [9] jquerylib_0.1.4   scales_1.3.0      yaml_2.3.10       fastmap_1.2.0    
-[13] lattice_0.22-6    ggplot2_3.5.1     R6_2.5.1          generics_0.1.3   
-[17] lmtest_0.9-40     MASS_7.3-61       tibble_3.2.1      munsell_0.5.1    
-[21] bslib_0.8.0       pillar_1.9.0      rlang_1.1.4       utf8_1.2.4       
-[25] cachem_1.1.0      xfun_0.49         sass_0.4.9        cli_3.6.3        
-[29] magrittr_2.0.3    digest_0.6.37     grid_4.4.2        mclust_6.1.1     
-[33] lifecycle_1.0.4   nlme_3.1-166      vctrs_0.6.5       evaluate_1.0.1   
-[37] glue_1.8.0        codetools_0.2-20  zoo_1.8-12        fansi_1.0.6      
-[41] colorspace_2.1-1  tools_4.4.2       pkgconfig_2.0.3   htmltools_0.5.8.1
+ [1] gtable_0.3.6      jsonlite_1.9.0    dplyr_1.1.4       compiler_4.5.0   
+ [5] tidyselect_1.2.1  gridExtra_2.3     jquerylib_0.1.4   systemfonts_1.2.1
+ [9] scales_1.3.0      textshaping_1.0.0 yaml_2.3.10       fastmap_1.2.0    
+[13] lattice_0.22-6    ggplot2_3.5.1     R6_2.6.1          generics_0.1.3   
+[17] lmtest_0.9-40     MASS_7.3-65       htmlwidgets_1.6.4 tibble_3.2.1     
+[21] desc_1.4.3        munsell_0.5.1     bslib_0.9.0       pillar_1.10.1    
+[25] rlang_1.1.5       cachem_1.1.0      xfun_0.51         fs_1.6.5         
+[29] sass_0.4.9        cli_3.6.4         pkgdown_2.1.1     magrittr_2.0.3   
+[33] digest_0.6.37     grid_4.5.0        mclust_6.1.1      lifecycle_1.0.4  
+[37] nlme_3.1-168      vctrs_0.6.5       evaluate_1.0.3    glue_1.8.0       
+[41] codetools_0.2-20  ragg_1.3.3        zoo_1.8-13        colorspace_2.1-1 
+[45] tools_4.5.0       pkgconfig_2.0.3   htmltools_0.5.8.1

Hardware info

Hardware info +

CPU model: AMD Ryzen 9 7950X 16-Core Processor

MemTotal:       64927788 kB

MemTotal:       64927780 kB

- - - - +

- - - + - - - - - + diff --git a/docs/articles/prebuilt/2022_dmta_pathway.html b/docs/articles/prebuilt/2022_dmta_pathway.html index c9f87cd8..d2ef55fb 100644 --- a/docs/articles/prebuilt/2022_dmta_pathway.html +++ b/docs/articles/prebuilt/2022_dmta_pathway.html @@ -1,396 +1,100 @@ - - - + - - - - - - - - - -Testing hierarchical pathway kinetics with residue data on dimethenamid and dimethenamid-P - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + +Testing hierarchical pathway kinetics with residue data on dimethenamid and dimethenamid-P • mkin + + + + + - + + + +

Testing hierarchical pathway kinetics with -residue data on dimethenamid and dimethenamid-P

Johannes Ranke

Last change on 20 April 2023, last compiled on 13 -Februar 2025

Introduction

+ + +

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, with parallel formation of two or more metabolites @@ -399,7 +103,7 @@ can be fitted with the mkin package.

173340 (Application of nonlinear hierarchical models to the kinetic evaluation of chemical degradation data) of the German Environment Agency carried out in 2022 and 2023.

The mkin package is used in version 1.2.9, which is currently under +

The mkin package is used in version 1.2.10, which is currently under development. It contains the test data, and the functions used in the evaluations. The saemix package is used as a backend for fitting the NLHM, but is also loaded to make the convergence plot @@ -408,26 +112,28 @@ function available.

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()
-
-# We need to start a new cluster after defining a compiled model that is
-# saved as a DLL to the user directory, therefore we define a function
-# This is used again after defining the pathway model
-start_cluster <- function(n_cores) {
-  if (Sys.info()["sysname"] == "Windows") {
-    ret <- makePSOCKcluster(n_cores)
-  } else {
-    ret <- makeForkCluster(n_cores)
-  }
-  return(ret)
-}

+library(mkin)
+library(knitr)
+library(saemix)
+library(parallel)
+n_cores <- detectCores()
+
+# We need to start a new cluster after defining a compiled model that is
+# saved as a DLL to the user directory, therefore we define a function
+# This is used again after defining the pathway model
+start_cluster <- function(n_cores) {
+  if (Sys.info()["sysname"] == "Windows") {
+    ret <- makePSOCKcluster(n_cores)
+  } else {
+    ret <- makeForkCluster(n_cores)
+  }
+  return(ret)
+}

Data

Data +

The test data are available in the mkin package as an object of class mkindsg (mkin dataset group) under the identifier dimethenamid_2018. The following preprocessing steps are @@ -448,41 +154,41 @@ and Elliot 2) are combined, resulting in dimethenamid (DMTA) data from six soils.

The following commented R code performs this preprocessing.

# Apply a function to each of the seven datasets in the mkindsg object to create a list
-dmta_ds <- lapply(1:7, function(i) {
-  ds_i <- dimethenamid_2018$ds[[i]]$data                     # Get a dataset
-  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"              # Rename DMTAP to DMTA
-  ds_i <- subset(ds_i, select = c("name", "time", "value")) # Select data
-  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]  # Normalise time
-  ds_i                                                       # Return the dataset
-})
-
-# Use dataset titles as names for the list elements
-names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-
-# Combine data for Elliot soil to obtain a named list with six elements
-dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) #
-dmta_ds[["Elliot 1"]] <- NULL
-dmta_ds[["Elliot 2"]] <- NULL

+# Apply a function to each of the seven datasets in the mkindsg object to create a list
+dmta_ds <- lapply(1:7, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data                     # Get a dataset
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"              # Rename DMTAP to DMTA
+  ds_i <- subset(ds_i, select = c("name", "time", "value")) # Select data
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]  # Normalise time
+  ds_i                                                       # Return the dataset
+})
+
+# Use dataset titles as names for the list elements
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+
+# Combine data for Elliot soil to obtain a named list with six elements
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) #
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL

The following tables show the 6 datasets.

for (ds_name in names(dmta_ds)) {
-  print(
-    kable(mkin_long_to_wide(dmta_ds[[ds_name]]),
-      caption = paste("Dataset", ds_name),
-      booktabs = TRUE, row.names = FALSE))
-    cat("\n\\clearpage\n")
-}

Hierarchical mkin fit of the SFO model with error model const
	npar	AIC	Chisq	Df	Pr(>Chisq)
f_saem_dfop_tc_no_ranef_k2

+for (ds_name in names(dmta_ds)) {
+  print(
+    kable(mkin_long_to_wide(dmta_ds[[ds_name]]),
+      caption = paste("Dataset", ds_name),
+      booktabs = TRUE, row.names = FALSE))
+    cat("\n\\clearpage\n")
+}

- - + - - + @@ -535,17 +241,15 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Calke
time	DMTA	M23	M27	M31
0

- +

- - + - - + @@ -661,17 +365,15 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Borstel
time	DMTA	M23	M27	M31
0.000000

- +

- - + - - + @@ -927,17 +629,15 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Flaach
time	DMTA	M23	M27	M31
0.0000000

- +

- - + - - + @@ -1109,17 +809,15 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset BBA 2.2
time	DMTA	M23	M27	M31
0.0000000

- +

- - + - - + @@ -1291,17 +989,15 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset BBA 2.3
time	DMTA	M23	M27	M31
0.0000000

- +

- - + - - + @@ -1642,8 +1338,9 @@ dmta_ds[["Elliot 2"]] <- NULL

Dataset Elliot
time	DMTA	M23	M27	M31
0.000000

Separate evaluations

Separate evaluations +

As a first step to obtain suitable starting parameters for the NLHM fits, we do separate fits of several variants of the pathway model used previously (Ranke et al. 2021), varying @@ -1651,71 +1348,71 @@ the kinetic model for the parent compound. Because the SFORB model often provides faster convergence than the DFOP model, and can sometimes be fitted where the DFOP model results in errors, it is included in the set of parent models tested here.

if (!dir.exists("dmta_dlls")) dir.create("dmta_dlls")
-m_sfo_path_1 <- mkinmod(
-  DMTA = mkinsub("SFO", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  name = "m_sfo_path", dll_dir = "dmta_dlls",
-  unload = TRUE, overwrite = TRUE,
-  quiet = TRUE
-)
-m_fomc_path_1 <- mkinmod(
-  DMTA = mkinsub("FOMC", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  name = "m_fomc_path", dll_dir = "dmta_dlls",
-  unload = TRUE, overwrite = TRUE,
-  quiet = TRUE
-)
-m_dfop_path_1 <- mkinmod(
-  DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  name = "m_dfop_path", dll_dir = "dmta_dlls",
-  unload = TRUE, overwrite = TRUE,
-  quiet = TRUE
-)
-m_sforb_path_1 <- mkinmod(
-  DMTA = mkinsub("SFORB", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  name = "m_sforb_path", dll_dir = "dmta_dlls",
-  unload = TRUE, overwrite = TRUE,
-  quiet = TRUE
-)
-m_hs_path_1 <- mkinmod(
-  DMTA = mkinsub("HS", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  name = "m_hs_path", dll_dir = "dmta_dlls",
-  unload = TRUE, overwrite = TRUE,
-  quiet = TRUE
-)
-cl <- start_cluster(n_cores)
-
-deg_mods_1 <- list(
-  sfo_path_1 = m_sfo_path_1,
-  fomc_path_1 = m_fomc_path_1,
-  dfop_path_1 = m_dfop_path_1,
-  sforb_path_1 = m_sforb_path_1,
-  hs_path_1 = m_hs_path_1)
-
-sep_1_const <- mmkin(
-  deg_mods_1,
-  dmta_ds,
-  error_model = "const",
-  quiet = TRUE)
-
-status(sep_1_const) |> kable()

- - - +

+if (!dir.exists("dmta_dlls")) dir.create("dmta_dlls")
+m_sfo_path_1 <- mkinmod(
+  DMTA = mkinsub("SFO", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  name = "m_sfo_path", dll_dir = "dmta_dlls",
+  unload = TRUE, overwrite = TRUE,
+  quiet = TRUE
+)
+m_fomc_path_1 <- mkinmod(
+  DMTA = mkinsub("FOMC", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  name = "m_fomc_path", dll_dir = "dmta_dlls",
+  unload = TRUE, overwrite = TRUE,
+  quiet = TRUE
+)
+m_dfop_path_1 <- mkinmod(
+  DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  name = "m_dfop_path", dll_dir = "dmta_dlls",
+  unload = TRUE, overwrite = TRUE,
+  quiet = TRUE
+)
+m_sforb_path_1 <- mkinmod(
+  DMTA = mkinsub("SFORB", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  name = "m_sforb_path", dll_dir = "dmta_dlls",
+  unload = TRUE, overwrite = TRUE,
+  quiet = TRUE
+)
+m_hs_path_1 <- mkinmod(
+  DMTA = mkinsub("HS", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  name = "m_hs_path", dll_dir = "dmta_dlls",
+  unload = TRUE, overwrite = TRUE,
+  quiet = TRUE
+)
+cl <- start_cluster(n_cores)
+
+deg_mods_1 <- list(
+  sfo_path_1 = m_sfo_path_1,
+  fomc_path_1 = m_fomc_path_1,
+  dfop_path_1 = m_dfop_path_1,
+  sforb_path_1 = m_sforb_path_1,
+  hs_path_1 = m_hs_path_1)
+
+sep_1_const <- mmkin(
+  deg_mods_1,
+  dmta_ds,
+  error_model = "const",
+  quiet = TRUE)
+
+status(sep_1_const) |> kable()

+ @@ -1723,8 +1420,7 @@ status(sep_1_const) |> kable() - - + @@ -1777,11 +1473,11 @@ status(sep_1_const) |> kable() constant variance converged (status OK). Most fits with DFOP or SFORB for the parent converged as well. The fits with HS for the parent did not converge with default settings.

sep_1_tc <- update(sep_1_const, error_model = "tc")
-status(sep_1_tc) |> kable()

	Calke	Borstel	BBA 2.2	BBA 2.3	Elliot
sfo_path_1

- - +

+sep_1_tc <- update(sep_1_const, error_model = "tc")
+status(sep_1_tc) |> kable()

+ @@ -1789,8 +1485,7 @@ status(sep_1_tc) |> kable() - - + @@ -1845,24 +1540,25 @@ different data sets when applying the DFOP and SFORB model and some additional convergence problems when using the FOMC model for the parent.

Hierarchichal model fits

Hierarchichal model fits +

The following code fits two sets of the corresponding hierarchical models to the data, one assuming constant variance, and one assuming two-component error.

saem_1 <- mhmkin(list(sep_1_const, sep_1_tc))

+saem_1 <- mhmkin(list(sep_1_const, sep_1_tc))

The run time for these fits was around two hours on five year old hardware. After a recent hardware upgrade these fits complete in less than twenty minutes.

status(saem_1) |> kable()

	Calke	Borstel	BBA 2.2	BBA 2.3	Elliot
sfo_path_1

- - +

+status(saem_1) |> kable()

+ - - + @@ -1893,17 +1589,16 @@ than twenty minutes.

	const	tc
sfo_path_1

According to the status function, all fits terminated successfully.

anova(saem_1) |> kable(digits = 1)

- - - +

+anova(saem_1) |> kable(digits = 1)

+ - - + @@ -1990,22 +1685,22 @@ assume a discontinuity, so the SFORB model is preferable from a mechanistic viewpoint. In addition, the information criteria AIC and BIC are very similar for HS and SFORB. Therefore, the SFORB model is selected here for further refinements.

Parameter identifiability based on the Fisher Information -Matrix

Parameter identifiability based on the Fisher Information +Matrix +

Using the illparms function, ill-defined statistical model parameters such as standard deviations of the degradation parameters in the population and error model parameters can be found.

illparms(saem_1) |> kable()

	npar	AIC	BIC	Lik
sfo_path_1 const

- - +

+illparms(saem_1) |> kable()

+ - - + @@ -2044,22 +1739,22 @@ two-component error, the random effect for the rate constant from reversibly bound DMTA to the free DMTA (k_DMTA_bound_free) is not well-defined. Therefore, the fit is updated without assuming a random effect for this parameter.

saem_sforb_path_1_tc_reduced <- update(saem_1[["sforb_path_1", "tc"]],
-  no_random_effect = "log_k_DMTA_bound_free")
-illparms(saem_sforb_path_1_tc_reduced)

+saem_sforb_path_1_tc_reduced <- update(saem_1[["sforb_path_1", "tc"]],
+  no_random_effect = "log_k_DMTA_bound_free")
+illparms(saem_sforb_path_1_tc_reduced)

As expected, no ill-defined parameters remain. The model comparison below shows that the reduced model is preferable.

anova(saem_1[["sforb_path_1", "tc"]], saem_sforb_path_1_tc_reduced) |> kable(digits = 1)

	const	tc
sfo_path_1

- - +

+anova(saem_1[["sforb_path_1", "tc"]], saem_sforb_path_1_tc_reduced) |> kable(digits = 1)

+ - - + @@ -2078,8 +1773,9 @@ below shows that the reduced model is preferable.

	npar	AIC	BIC	Lik
saem_sforb_path_1_tc_reduced

The convergence plot of the refined fit is shown below.

plot(saem_sforb_path_1_tc_reduced$so, plot.type = "convergence")

+plot(saem_sforb_path_1_tc_reduced$so, plot.type = "convergence")

For some parameters, for example for f_DMTA_ilr_1 and f_DMTA_ilr_2, i.e. for two of the parameters determining the formation fractions of the parallel formation of the three @@ -2091,15 +1787,18 @@ the parameter estimates very much, and it is proposed that the fit is acceptable. No numeric convergence criterion is implemented in saemix.

Alternative check of parameter identifiability

Alternative check of parameter identifiability +

As an alternative check of parameter identifiability (Duchesne et al. 2021), multistart runs were performed on the basis of the refined fit shown above.

saem_sforb_path_1_tc_reduced_multi <- multistart(saem_sforb_path_1_tc_reduced,
-  n = 32, cores = 10)

+saem_sforb_path_1_tc_reduced_multi <- multistart(saem_sforb_path_1_tc_reduced,
+  n = 32, cores = 10)

 
   (subscript) logical subscript too long

print(saem_sforb_path_1_tc_reduced_multi)

+print(saem_sforb_path_1_tc_reduced_multi)

<multistart> object with 32 fits:
  E OK 
  7 25 
@@ -2112,11 +1811,11 @@ the SAEM algorithm leads to parameter combinations for the degradation
 model that the numerical integration routine cannot cope with. Because
 of this variation of initial parameters, some of the model fits take up
 to two times more time than the original fit.
-par(mar = c(12.1, 4.1, 2.1, 2.1))
-parplot(saem_sforb_path_1_tc_reduced_multi, ylim = c(0.5, 2), las = 2)
++par(mar = c(12.1, 4.1, 2.1, 2.1))
+parplot(saem_sforb_path_1_tc_reduced_multi, ylim = c(0.5, 2), las = 2)
 
-
-
+

 Parameter boxplots for the multistart runs that succeeded
 
 
@@ -2126,103 +1825,114 @@ independent of the starting parameters, and there are no remaining
 ill-defined parameters.


-
-Plots of selected fits
+
+Plots of selected fits
+
 The SFORB pathway fits with full and reduced parameter distribution
 model are shown below.
-plot(saem_1[["sforb_path_1", "tc"]])
++plot(saem_1[["sforb_path_1", "tc"]])
 
-
-
+

 SFORB pathway fit with two-component error
 
 
-plot(saem_sforb_path_1_tc_reduced)
++plot(saem_sforb_path_1_tc_reduced)
 
-
-
+

 SFORB pathway fit with two-component error, reduced parameter model
 
 
 Plots of the remaining fits and listings for all successful fits are
 shown in the Appendix.
-stopCluster(cl)
++stopCluster(cl)
 
-
-Conclusions
+
+Conclusions
+
 Pathway fits with SFO, FOMC, DFOP, SFORB and HS models for the parent
 compound could be successfully performed.
 
-
-Acknowledgements
+
+Acknowledgements
+
 The helpful comments by Janina Wöltjen of the German Environment
 Agency on earlier versions of this document are gratefully
 acknowledged.
 
-
-References
+
+References
+
 
 
 Duchesne, Ronan, Anissa Guillemin, Olivier Gandrillon, and Fabien
 Crauste. 2021. “Practical Identifiability in the Frame of
 Nonlinear Mixed Effects Models: The Example of the in Vitro
-Erythropoiesis.” BMC Bioinformatics 22 (478). https://doi.org/10.1186/s12859-021-04373-4.
+Erythropoiesis.” BMC Bioinformatics 22 (478). https://doi.org/10.1186/s12859-021-04373-4.
 
 
 Ranke, Johannes, Janina Wöltjen, Jana Schmidt, and Emmanuelle Comets.
 2021. “Taking Kinetic Evaluations of Degradation Data to the Next
 Level with Nonlinear Mixed-Effects Models.” Environments
-8 (8). https://doi.org/10.3390/environments8080071.
+8 (8). https://doi.org/10.3390/environments8080071.
 
 
 
-
-Appendix
-
-Plots of hierarchical fits not selected for refinement
-plot(saem_1[["sfo_path_1", "tc"]])
+
+Appendix
+
+
+Plots of hierarchical fits not selected for refinement
+
++plot(saem_1[["sfo_path_1", "tc"]])
 
-
-
+

 SFO pathway fit with two-component error
 
 
-plot(saem_1[["fomc_path_1", "tc"]])
++plot(saem_1[["fomc_path_1", "tc"]])
 
-
-
+

 FOMC pathway fit with two-component error
 
 
-plot(saem_1[["sforb_path_1", "tc"]])
++plot(saem_1[["sforb_path_1", "tc"]])
 
-
-
+

 HS pathway fit with two-component error
 
 
 
-
-Hierarchical model fit listings
-
-Fits with random effects for all degradation parameters
+
+Hierarchical model fit listings
+
+
+Fits with random effects for all degradation parameters
+
 
 
-
-Improved fit of the SFORB pathway model with two-component
-error
+
+Improved fit of the SFORB pathway model with two-component
+error
+
 
 
 
-
-Session info
-R version 4.4.2 (2024-10-31)
+
+Session info
+
+R version 4.5.0 (2025-04-11)
 Platform: x86_64-pc-linux-gnu
 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
+LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0  LAPACK version 3.11.0
 
 locale:
  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
@@ -2240,75 +1950,52 @@ attached base packages:
 [8] base     
 
 other attached packages:
-[1] saemix_3.3      npde_3.5        knitr_1.49      mkin_1.2.9     
-[5] rmarkdown_2.29  nvimcom_0.9-167
+[1] rmarkdown_2.29  nvimcom_0.9-167 saemix_3.3      npde_3.5       
+[5] knitr_1.49      mkin_1.2.10    
 
 loaded via a namespace (and not attached):
- [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    lattice_0.22-6   
- [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.1    grid_4.4.2       
- [9] fastmap_1.2.0     jsonlite_1.8.9    processx_3.8.4    pkgbuild_1.4.5   
-[13] deSolve_1.40      mclust_6.1.1      ps_1.8.1          gridExtra_2.3    
-[17] fansi_1.0.6       scales_1.3.0      codetools_0.2-20  jquerylib_0.1.4  
-[21] cli_3.6.3         rlang_1.1.4       munsell_0.5.1     cachem_1.1.0     
-[25] yaml_2.3.10       tools_4.4.2       inline_0.3.20     colorout_1.3-2   
-[29] dplyr_1.1.4       colorspace_2.1-1  ggplot2_3.5.1     vctrs_0.6.5      
-[33] R6_2.5.1          zoo_1.8-12        lifecycle_1.0.4   MASS_7.3-61      
-[37] pkgconfig_2.0.3   callr_3.7.6       pillar_1.9.0      bslib_0.8.0      
-[41] gtable_0.3.6      glue_1.8.0        xfun_0.49         tibble_3.2.1     
-[45] lmtest_0.9-40     tidyselect_1.2.1  htmltools_0.5.8.1 nlme_3.1-166     
-[49] compiler_4.4.2   
+ [1] sass_0.4.9        generics_0.1.3    lattice_0.22-6    digest_0.6.37    
+ [5] magrittr_2.0.3    evaluate_1.0.3    grid_4.5.0        fastmap_1.2.0    
+ [9] jsonlite_1.9.0    processx_3.8.6    pkgbuild_1.4.6    deSolve_1.40     
+[13] mclust_6.1.1      ps_1.9.0          gridExtra_2.3     scales_1.3.0     
+[17] codetools_0.2-20  textshaping_1.0.0 jquerylib_0.1.4   cli_3.6.4        
+[21] rlang_1.1.5       munsell_0.5.1     cachem_1.1.0      yaml_2.3.10      
+[25] tools_4.5.0       inline_0.3.21     dplyr_1.1.4       colorspace_2.1-1 
+[29] ggplot2_3.5.1     vctrs_0.6.5       R6_2.6.1          zoo_1.8-13       
+[33] lifecycle_1.0.4   fs_1.6.5          htmlwidgets_1.6.4 MASS_7.3-65      
+[37] ragg_1.3.3        pkgconfig_2.0.3   desc_1.4.3        callr_3.7.6      
+[41] pkgdown_2.1.1     pillar_1.10.1     bslib_0.9.0       gtable_0.3.6     
+[45] glue_1.8.0        systemfonts_1.2.1 xfun_0.51         tibble_3.2.1     
+[49] lmtest_0.9-40     tidyselect_1.2.1  htmltools_0.5.8.1 nlme_3.1-168     
+[53] compiler_4.5.0   
 
-
-Hardware info
+
+Hardware info
+
 CPU model: AMD Ryzen 9 7950X 16-Core Processor
-MemTotal:       64927788 kB
+MemTotal:       64927780 kB
 
 
-
-
-
-
+

1.949	2.411
Ryzen 9 7950X	Linux	1.2.10	3.3	1.115	2.263	1.991	2.193
Ryzen 9 7950X	Linux	1.2.10	3.3	1.876	1.922	2.423	1.849

3.256	3.322
Ryzen 9 7950X	Linux	1.2.10	3.3	2.380	3.147	3.131	3.331
Ryzen 9 7950X	Linux	1.2.10	3.3	2.397	3.443	3.414	3.311

12.420	289.338
Ryzen 9 7950X	Linux	1.2.10	3.3	11.581	296.184
Ryzen 9 7950X	Linux	1.2.10	3.3	12.067	278.742

CPU	3.3	485.836
Ryzen 9 7950X	Linux	1.2.10	3.3	486.856
Ryzen 9 7950X	Linux	1.2.10	3.3	474.550

Last change 31 January 2019 -(rebuilt 2025-02-13)

Last change 18 May 2023 -(rebuilt 2025-02-13)

Example evaluations with (generalised) nonlinear least squares

Miscellaneous

Last change 18 May 2023 -(rebuilt 2025-02-13)

Abstract

Abstract +

Background

Background +

Derived software tools

Derived software tools +

Unique features

Unique features +

Internal parameter transformations

Internal parameter transformations +

Confidence intervals based on transformed parameters

Confidence intervals based on transformed parameters +

Parameter t-test based on untransformed parameters

Parameter t-test based on untransformed parameters +

References

References +

Testing hierarchical pathway kinetics with -residue data on cyantraniliprole

Johannes Ranke

Last change on 13 February 2023, last compiled on 13 -Februar 2025

Testing hierarchical pathway kinetics with residue data on cyantraniliprole

Johannes +Ranke

Last change on 13 February +2023, last compiled on 12 Mai 2025

Introduction

Introduction +

Test data

Test data +

Parent only evaluations

Parent only evaluations +

Pathway fits

Evaluations with pathway established previously

Pathway fits +

Evaluations with pathway established previously +

Alternative pathway fits

Alternative pathway fits +

Refinement of alternative pathway fits

Refinement of alternative pathway fits +

Conclusion

Conclusion +

Acknowledgements

Acknowledgements +

Appendix

Plots of fits that were not refined further

Appendix +

Plots of fits that were not refined further +

Hierarchical fit listings

Pathway 1

Hierarchical fit listings +

Pathway 1 +

Pathway 2

Pathway 2 +

Pathway 2, refined fits

Pathway 2, refined fits +

Session info

Session info +

Hardware info

Hardware info +

Separate evaluations

Separate evaluations +

Hierarchichal model fits

Hierarchichal model fits +

Parameter identifiability based on the Fisher Information -Matrix

Parameter identifiability based on the Fisher Information +Matrix +

Alternative check of parameter identifiability

Alternative check of parameter identifiability +

Conclusions

Conclusions +

Acknowledgements

Acknowledgements +

References

References +

Appendix

Hierarchical model fit listings

Appendix +

Hierarchical model fit listings +

Hierarchical model convergence plots