From 91a5834dd701211f929fd25419dc34561ce3b4e7 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 14 Feb 2025 09:15:20 +0100 Subject: Initialize dev docs --- docs/dev/articles/web_only/FOCUS_Z.html | 446 ++++++++ .../figure-html/FOCUS_2006_Z_fits_1-1.png | Bin 0 -> 67444 bytes .../figure-html/FOCUS_2006_Z_fits_10-1.png | Bin 0 -> 107009 bytes .../figure-html/FOCUS_2006_Z_fits_11-1.png | Bin 0 -> 106098 bytes .../figure-html/FOCUS_2006_Z_fits_11a-1.png | Bin 0 -> 76831 bytes .../figure-html/FOCUS_2006_Z_fits_11b-1.png | Bin 0 -> 37494 bytes .../figure-html/FOCUS_2006_Z_fits_2-1.png | Bin 0 -> 67444 bytes .../figure-html/FOCUS_2006_Z_fits_3-1.png | Bin 0 -> 67259 bytes .../figure-html/FOCUS_2006_Z_fits_5-1.png | Bin 0 -> 81751 bytes .../figure-html/FOCUS_2006_Z_fits_6-1.png | Bin 0 -> 106406 bytes .../figure-html/FOCUS_2006_Z_fits_7-1.png | Bin 0 -> 105359 bytes .../figure-html/FOCUS_2006_Z_fits_9-1.png | Bin 0 -> 89110 bytes docs/dev/articles/web_only/NAFTA_examples.html | 1049 ++++++++++++++++++ .../NAFTA_examples_files/figure-html/p10-1.png | Bin 0 -> 81542 bytes .../NAFTA_examples_files/figure-html/p11-1.png | Bin 0 -> 77600 bytes .../NAFTA_examples_files/figure-html/p12a-1.png | Bin 0 -> 83512 bytes .../NAFTA_examples_files/figure-html/p12b-1.png | Bin 0 -> 71409 bytes .../NAFTA_examples_files/figure-html/p13-1.png | Bin 0 -> 78613 bytes .../NAFTA_examples_files/figure-html/p14-1.png | Bin 0 -> 81557 bytes .../NAFTA_examples_files/figure-html/p15a-1.png | Bin 0 -> 77580 bytes .../NAFTA_examples_files/figure-html/p15b-1.png | Bin 0 -> 79748 bytes .../NAFTA_examples_files/figure-html/p16-1.png | Bin 0 -> 94540 bytes .../NAFTA_examples_files/figure-html/p5a-1.png | Bin 0 -> 83061 bytes .../NAFTA_examples_files/figure-html/p5b-1.png | Bin 0 -> 81186 bytes .../NAFTA_examples_files/figure-html/p6-1.png | Bin 0 -> 83142 bytes .../NAFTA_examples_files/figure-html/p7-1.png | Bin 0 -> 102934 bytes .../NAFTA_examples_files/figure-html/p8-1.png | Bin 0 -> 93632 bytes .../NAFTA_examples_files/figure-html/p9a-1.png | Bin 0 -> 78782 bytes .../NAFTA_examples_files/figure-html/p9b-1.png | Bin 0 -> 76883 bytes docs/dev/articles/web_only/benchmarks.html | 1133 ++++++++++++++++++++ docs/dev/articles/web_only/compiled_models.html | 231 ++++ docs/dev/articles/web_only/dimethenamid_2018.html | 651 +++++++++++ .../figure-html/f_parent_mkin_dfop_const-1.png | Bin 0 -> 57786 bytes .../f_parent_mkin_dfop_const_test-1.png | Bin 0 -> 57786 bytes .../figure-html/f_parent_mkin_dfop_tc_test-1.png | Bin 0 -> 59146 bytes .../figure-html/f_parent_mkin_sfo_const-1.png | Bin 0 -> 55982 bytes .../figure-html/f_parent_saemix_dfop_const-1.png | Bin 0 -> 36337 bytes .../figure-html/f_parent_saemix_dfop_tc-1.png | Bin 0 -> 29343 bytes .../figure-html/f_parent_saemix_sfo_const-1.png | Bin 0 -> 32949 bytes .../figure-html/f_parent_saemix_sfo_tc-1.png | Bin 0 -> 28538 bytes .../figure-html/plot_parent_nlme-1.png | Bin 0 -> 59209 bytes docs/dev/articles/web_only/multistart.html | 179 ++++ .../figure-html/unnamed-chunk-3-1.png | Bin 0 -> 62571 bytes .../figure-html/unnamed-chunk-4-1.png | Bin 0 -> 53560 bytes .../figure-html/unnamed-chunk-5-1.png | Bin 0 -> 22804 bytes docs/dev/articles/web_only/saem_benchmarks.html | 789 ++++++++++++++ 46 files changed, 4478 insertions(+) create mode 100644 docs/dev/articles/web_only/FOCUS_Z.html create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_1-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_2-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_3-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_5-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_7-1.png create mode 100644 docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples.html create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p10-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p11-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12a-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12b-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p13-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15a-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15b-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p16-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5a-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5b-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p6-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p7-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p8-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9a-1.png create mode 100644 docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png create mode 100644 docs/dev/articles/web_only/benchmarks.html create mode 100644 docs/dev/articles/web_only/compiled_models.html create mode 100644 docs/dev/articles/web_only/dimethenamid_2018.html create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const_test-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_tc_test-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_sfo_const-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_const-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_tc-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_const-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_tc-1.png create mode 100644 docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/plot_parent_nlme-1.png create mode 100644 docs/dev/articles/web_only/multistart.html create mode 100644 docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png create mode 100644 docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-4-1.png create mode 100644 docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-5-1.png create mode 100644 docs/dev/articles/web_only/saem_benchmarks.html (limited to 'docs/dev/articles/web_only') diff --git a/docs/dev/articles/web_only/FOCUS_Z.html b/docs/dev/articles/web_only/FOCUS_Z.html new file mode 100644 index 00000000..40ce9e01 --- /dev/null +++ b/docs/dev/articles/web_only/FOCUS_Z.html @@ -0,0 +1,446 @@ + + + + + + + +Example evaluation of FOCUS dataset Z • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Example evaluation of FOCUS dataset Z

+

Johannes +Ranke

+ +

Last change 16 January 2018 +(rebuilt 2025-02-14)

+ + Source: vignettes/web_only/FOCUS_Z.rmd +
FOCUS_Z.rmd
+
+ + + +

Wissenschaftlicher Berater, Kronacher +Str. 12, 79639 Grenzach-Wyhlen, Germany
Privatdozent at the +University of Bremen

+
+

The data +

+

The following code defines the example dataset from Appendix 7 to the +FOCUS kinetics report (FOCUS Work Group on +Degradation Kinetics 2014, 354).

+
+library(mkin, quietly = TRUE)
+LOD = 0.5
+FOCUS_2006_Z = data.frame(
+  t = c(0, 0.04, 0.125, 0.29, 0.54, 1, 2, 3, 4, 7, 10, 14, 21,
+        42, 61, 96, 124),
+  Z0 = c(100, 81.7, 70.4, 51.1, 41.2, 6.6, 4.6, 3.9, 4.6, 4.3, 6.8,
+         2.9, 3.5, 5.3, 4.4, 1.2, 0.7),
+  Z1 = c(0, 18.3, 29.6, 46.3, 55.1, 65.7, 39.1, 36, 15.3, 5.6, 1.1,
+         1.6, 0.6, 0.5 * LOD, NA, NA, NA),
+  Z2 = c(0, NA, 0.5 * LOD, 2.6, 3.8, 15.3, 37.2, 31.7, 35.6, 14.5,
+         0.8, 2.1, 1.9, 0.5 * LOD, NA, NA, NA),
+  Z3 = c(0, NA, NA, NA, NA, 0.5 * LOD, 9.2, 13.1, 22.3, 28.4, 32.5,
+         25.2, 17.2, 4.8, 4.5, 2.8, 4.4))
+
+FOCUS_2006_Z_mkin <- mkin_wide_to_long(FOCUS_2006_Z)
+
+
+

Parent and one metabolite +

+

The next step is to set up the models used for the kinetic analysis. +As the simultaneous fit of parent and the first metabolite is usually +straightforward, Step 1 (SFO for parent only) is skipped here. We start +with the model 2a, with formation and decline of metabolite Z1 and the +pathway from parent directly to sink included (default in mkin).

+
+Z.2a <- mkinmod(Z0 = mkinsub("SFO", "Z1"),
+                Z1 = mkinsub("SFO"))
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.2a <- mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE): Observations with
+## value of zero were removed from the data
+
+plot_sep(m.Z.2a)
+

+
+summary(m.Z.2a, data = FALSE)$bpar
+
##            Estimate se_notrans t value     Pr(>t)    Lower    Upper
+## Z0_0       97.01488   3.301084 29.3888 3.2971e-21 91.66556 102.3642
+## k_Z0        2.23601   0.207078 10.7979 3.3309e-11  1.95303   2.5600
+## k_Z1        0.48212   0.063265  7.6207 2.8154e-08  0.40341   0.5762
+## f_Z0_to_Z1  1.00000   0.094764 10.5525 5.3560e-11  0.00000   1.0000
+## sigma       4.80411   0.635638  7.5579 3.2592e-08  3.52677   6.0815
+

As obvious from the parameter summary (the component of the summary), +the kinetic rate constant from parent compound Z to sink is very small +and the t-test for this parameter suggests that it is not significantly +different from zero. This suggests, in agreement with the analysis in +the FOCUS kinetics report, to simplify the model by removing the pathway +to sink.

+

A similar result can be obtained when formation fractions are used in +the model formulation:

+
+Z.2a.ff <- mkinmod(Z0 = mkinsub("SFO", "Z1"),
+                   Z1 = mkinsub("SFO"),
+                   use_of_ff = "max")
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.2a.ff <- mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE): Observations with
+## value of zero were removed from the data
+
+plot_sep(m.Z.2a.ff)
+

+
+summary(m.Z.2a.ff, data = FALSE)$bpar
+
##            Estimate se_notrans t value     Pr(>t)    Lower    Upper
+## Z0_0       97.01488   3.301084 29.3888 3.2971e-21 91.66556 102.3642
+## k_Z0        2.23601   0.207078 10.7979 3.3309e-11  1.95303   2.5600
+## k_Z1        0.48212   0.063265  7.6207 2.8154e-08  0.40341   0.5762
+## f_Z0_to_Z1  1.00000   0.094764 10.5525 5.3560e-11  0.00000   1.0000
+## sigma       4.80411   0.635638  7.5579 3.2592e-08  3.52677   6.0815
+

Here, the ilr transformed formation fraction fitted in the model +takes a very large value, and the backtransformed formation fraction +from parent Z to Z1 is practically unity. Here, the covariance matrix +used for the calculation of confidence intervals is not returned as the +model is overparameterised.

+

A simplified model is obtained by removing the pathway to the sink. +

+

In the following, we use the parameterisation with formation +fractions in order to be able to compare with the results in the FOCUS +guidance, and as it makes it easier to use parameters obtained in a +previous fit when adding a further metabolite.

+
+Z.3 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
+               Z1 = mkinsub("SFO"), use_of_ff = "max")
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE): Observations with
+## value of zero were removed from the data
+
+plot_sep(m.Z.3)
+

+
+summary(m.Z.3, data = FALSE)$bpar
+
##       Estimate se_notrans t value     Pr(>t)    Lower    Upper
+## Z0_0  97.01488   2.597342  37.352 2.0106e-24 91.67597 102.3538
+## k_Z0   2.23601   0.146904  15.221 9.1477e-15  1.95354   2.5593
+## k_Z1   0.48212   0.041727  11.554 4.8268e-12  0.40355   0.5760
+## sigma  4.80411   0.620208   7.746 1.6110e-08  3.52925   6.0790
+

As there is only one transformation product for Z0 and no pathway to +sink, the formation fraction is internally fixed to unity.

+
+
+

Metabolites Z2 and Z3 +

+

As suggested in the FOCUS report, the pathway to sink was removed for +metabolite Z1 as well in the next step. While this step appears +questionable on the basis of the above results, it is followed here for +the purpose of comparison. Also, in the FOCUS report, it is assumed that +there is additional empirical evidence that Z1 quickly and exclusively +hydrolyses to Z2.

+
+Z.5 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
+               Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+               Z2 = mkinsub("SFO"), use_of_ff = "max")
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.5 <- mkinfit(Z.5, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.5, FOCUS_2006_Z_mkin, quiet = TRUE): Observations with
+## value of zero were removed from the data
+
+plot_sep(m.Z.5)
+

+

Finally, metabolite Z3 is added to the model. We use the optimised +differential equation parameter values from the previous fit in order to +accelerate the optimization.

+
+Z.FOCUS <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
+                   Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                   Z2 = mkinsub("SFO", "Z3"),
+                   Z3 = mkinsub("SFO"),
+                   use_of_ff = "max")
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin,
+                     parms.ini = m.Z.5$bparms.ode,
+                     quiet = TRUE)
+
## Warning in mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = m.Z.5$bparms.ode, :
+## Observations with value of zero were removed from the data
+
## Warning in mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = m.Z.5$bparms.ode, : Optimisation did not converge:
+## false convergence (8)
+
+plot_sep(m.Z.FOCUS)
+

+
+summary(m.Z.FOCUS, data = FALSE)$bpar
+
##             Estimate se_notrans t value     Pr(>t)     Lower      Upper
+## Z0_0       96.842440   1.994291 48.5598 4.0226e-42 92.830421 100.854459
+## k_Z0        2.215425   0.118457 18.7023 1.0404e-23  1.989490   2.467019
+## k_Z1        0.478307   0.028257 16.9272 6.2332e-22  0.424709   0.538669
+## k_Z2        0.451642   0.042139 10.7178 1.6304e-14  0.374348   0.544894
+## k_Z3        0.058692   0.015245  3.8499 1.7803e-04  0.034804   0.098975
+## f_Z2_to_Z3  0.471483   0.058348  8.0806 9.6585e-11  0.357720   0.588287
+## sigma       3.984431   0.383402 10.3923 4.5576e-14  3.213126   4.755737
+
+endpoints(m.Z.FOCUS)
+
## $ff
+##   Z2_Z3 Z2_sink 
+## 0.47148 0.52852 
+## 
+## $distimes
+##        DT50    DT90
+## Z0  0.31287  1.0393
+## Z1  1.44917  4.8140
+## Z2  1.53473  5.0983
+## Z3 11.80991 39.2317
+

This fit corresponds to the final result chosen in Appendix 7 of the +FOCUS report. Confidence intervals returned by mkin are based on +internally transformed parameters, however.

+
+
+

Using the SFORB model +

+

As the FOCUS report states, there is a certain tailing of the time +course of metabolite Z3. Also, the time course of the parent compound is +not fitted very well using the SFO model, as residues at a certain low +level remain.

+

Therefore, an additional model is offered here, using the single +first-order reversible binding (SFORB) model for metabolite Z3. As +expected, the +χ2\chi^2 +error level is lower for metabolite Z3 using this model and the +graphical fit for Z3 is improved. However, the covariance matrix is not +returned.

+
+Z.mkin.1 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFORB"))
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE): Observations
+## with value of zero were removed from the data
+
+plot_sep(m.Z.mkin.1)
+

+
+summary(m.Z.mkin.1, data = FALSE)$cov.unscaled
+
## NULL
+

Therefore, a further stepwise model building is performed starting +from the stage of parent and two metabolites, starting from the +assumption that the model fit for the parent compound can be improved by +using the SFORB model.

+
+Z.mkin.3 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO"))
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE)
+
## Warning in mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE): Observations
+## with value of zero were removed from the data
+
+plot_sep(m.Z.mkin.3)
+

+

This results in a much better representation of the behaviour of the +parent compound Z0.

+

Finally, Z3 is added as well. These models appear overparameterised +(no covariance matrix returned) if the sink for Z1 is left in the +models.

+
+Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFO"))
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin,
+                      parms.ini = m.Z.mkin.3$bparms.ode,
+                      quiet = TRUE)
+
## Warning in mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, parms.ini =
+## m.Z.mkin.3$bparms.ode, : Observations with value of zero were removed from the
+## data
+
+plot_sep(m.Z.mkin.4)
+

+

The error level of the fit, but especially of metabolite Z3, can be +improved if the SFORB model is chosen for this metabolite, as this model +is capable of representing the tailing of the metabolite decline +phase.

+
+Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFORB"))
+
## Temporary DLL for differentials generated and loaded
+
+m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin,
+                      parms.ini = m.Z.mkin.4$bparms.ode[1:4],
+                      quiet = TRUE)
+
## Warning in mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini =
+## m.Z.mkin.4$bparms.ode[1:4], : Observations with value of zero were removed from
+## the data
+
+plot_sep(m.Z.mkin.5)
+

+

The summary view of the backtransformed parameters shows that we get +no confidence intervals due to overparameterisation. As the optimized is +excessively small, it seems reasonable to fix it to zero.

+
+m.Z.mkin.5a <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin,
+                       parms.ini = c(m.Z.mkin.5$bparms.ode[1:7],
+                                     k_Z3_bound_free = 0),
+                       fixed_parms = "k_Z3_bound_free",
+                       quiet = TRUE)
+
## Warning in mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini =
+## c(m.Z.mkin.5$bparms.ode[1:7], : Observations with value of zero were removed
+## from the data
+
+plot_sep(m.Z.mkin.5a)
+

+

As expected, the residual plots for Z0 and Z3 are more random than in +the case of the all SFO model for which they were shown above. In +conclusion, the model is proposed as the best-fit model for the dataset +from Appendix 7 of the FOCUS report.

+

A graphical representation of the confidence intervals can finally be +obtained.

+
+mkinparplot(m.Z.mkin.5a)
+

+

The endpoints obtained with this model are

+
+endpoints(m.Z.mkin.5a)
+
## $ff
+## Z0_free   Z2_Z3 Z2_sink Z3_free 
+## 1.00000 0.53656 0.46344 1.00000 
+## 
+## $SFORB
+##     Z0_b1     Z0_b2      Z0_g     Z3_b1     Z3_b2      Z3_g 
+## 2.4471342 0.0075124 0.9519866 0.0800071 0.0000000 0.9347816 
+## 
+## $distimes
+##      DT50   DT90 DT50back DT50_Z0_b1 DT50_Z0_b2 DT50_Z3_b1 DT50_Z3_b2
+## Z0 0.3043 1.1848  0.35666    0.28325     92.267         NA         NA
+## Z1 1.5148 5.0320       NA         NA         NA         NA         NA
+## Z2 1.6414 5.4526       NA         NA         NA         NA         NA
+## Z3     NA     NA       NA         NA         NA     8.6636        Inf
+

It is clear the degradation rate of Z3 towards the end of the +experiment is very low as DT50_Z3_b2 (the second Eigenvalue of the +system of two differential equations representing the SFORB system for +Z3, corresponding to the slower rate constant of the DFOP model) is +reported to be infinity. However, this appears to be a feature of the +data.

+
+
+

References +

+ +
+
+FOCUS Work Group on Degradation Kinetics. 2014. Generic Guidance for +Estimating Persistence and Degradation Kinetics from Environmental Fate +Studies on Pesticides in EU Registration. 1.1 ed. http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics. +
+
+
+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_1-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_1-1.png new file mode 100644 index 00000000..98bc135b Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_1-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png new file mode 100644 index 00000000..c1011a35 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png new file mode 100644 index 00000000..dfd2dd50 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png new file mode 100644 index 00000000..74173f36 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png new file mode 100644 index 00000000..1c5793cc Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_2-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_2-1.png new file mode 100644 index 00000000..98bc135b Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_2-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_3-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_3-1.png new file mode 100644 index 00000000..0380ba43 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_3-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_5-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_5-1.png new file mode 100644 index 00000000..8c594ec9 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_5-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png new file mode 100644 index 00000000..84d473d6 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_7-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_7-1.png new file mode 100644 index 00000000..87af8874 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_7-1.png differ diff --git a/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png new file mode 100644 index 00000000..492cdcc8 Binary files /dev/null and b/docs/dev/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples.html b/docs/dev/articles/web_only/NAFTA_examples.html new file mode 100644 index 00000000..6bf39b2c --- /dev/null +++ b/docs/dev/articles/web_only/NAFTA_examples.html @@ -0,0 +1,1049 @@ + + + + + + + +Evaluation of example datasets from Attachment 1 to the US EPA SOP for the NAFTA guidance • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Evaluation of example datasets from Attachment 1 to the US EPA SOP for the NAFTA guidance

+

Johannes +Ranke

+ +

26 February 2019 (rebuilt +2025-02-14)

+ + Source: vignettes/web_only/NAFTA_examples.rmd +
NAFTA_examples.rmd
+
+ + + +
+

Introduction +

+

In this document, the example evaluations provided in Attachment 1 to +the SOP of US EPA for using the NAFTA guidance (US EPA 2015) are repeated using mkin. The +original evaluations reported in the attachment were performed using +PestDF in version 0.8.4. Note that PestDF 0.8.13 is the version +distributed at the US EPA website today (2019-02-26).

+

The datasets are now distributed with the mkin package.

+
+
+

Examples where DFOP did not converge with PestDF 0.8.4 +

+

In attachment 1, it is reported that the DFOP model does not converge +for these datasets when PestDF 0.8.4 was used. For all four datasets, +the DFOP model can be fitted with mkin (see below). The negative +half-life given by PestDF 0.8.4 for these fits appears to be the result +of a bug. The results for the other two models (SFO and IORE) are the +same.

+
+

Example on page 5, upper panel +

+
+p5a <- nafta(NAFTA_SOP_Attachment[["p5a"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p5a)
+

+
+print(p5a)
+
## Sums of squares:
+##       SFO      IORE      DFOP 
+## 465.21753  56.27506  32.06401 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 64.4304
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)  Lower   Upper
+## parent_0  95.8401 4.67e-21 92.245 99.4357
+## k_parent   0.0102 3.92e-12  0.009  0.0117
+## sigma      4.8230 3.81e-06  3.214  6.4318
+## 
+## $IORE
+##                Estimate Pr(>t)    Lower    Upper
+## parent_0       1.01e+02     NA 9.91e+01 1.02e+02
+## k__iore_parent 1.54e-05     NA 4.08e-06 5.84e-05
+## N_parent       2.57e+00     NA 2.25e+00 2.89e+00
+## sigma          1.68e+00     NA 1.12e+00 2.24e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower    Upper
+## parent_0 9.99e+01 1.41e-26 98.8116 101.0810
+## k1       2.67e-02 5.05e-06  0.0243   0.0295
+## k2       3.41e-12 5.00e-01  0.0000      Inf
+## g        6.47e-01 3.67e-06  0.6248   0.6677
+## sigma    1.27e+00 8.91e-06  0.8395   1.6929
+## 
+## 
+## DTx values:
+##      DT50     DT90 DT50_rep
+## SFO  67.7 2.25e+02 6.77e+01
+## IORE 58.2 1.07e+03 3.22e+02
+## DFOP 55.5 3.70e+11 2.03e+11
+## 
+## Representative half-life:
+## [1] 321.51
+
+
+

Example on page 5, lower panel +

+
+p5b <- nafta(NAFTA_SOP_Attachment[["p5b"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p5b)
+

+
+print(p5b)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 94.81123 10.10936  7.55871 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 11.77879
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0   96.497 2.32e-24 94.85271 98.14155
+## k_parent    0.008 3.42e-14  0.00737  0.00869
+## sigma       2.295 1.22e-05  1.47976  3.11036
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower    Upper
+## parent_0       9.85e+01 1.17e-28 9.79e+01 9.92e+01
+## k__iore_parent 1.53e-04 6.50e-03 7.21e-05 3.26e-04
+## N_parent       1.94e+00 5.88e-13 1.76e+00 2.12e+00
+## sigma          7.49e-01 1.63e-05 4.82e-01 1.02e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0 9.84e+01 1.24e-27 97.8078 98.9187
+## k1       1.55e-02 4.10e-04  0.0143  0.0167
+## k2       9.07e-12 5.00e-01  0.0000     Inf
+## g        6.89e-01 2.92e-03  0.6626  0.7142
+## sigma    6.48e-01 2.38e-05  0.4147  0.8813
+## 
+## 
+## DTx values:
+##      DT50     DT90 DT50_rep
+## SFO  86.6 2.88e+02 8.66e+01
+## IORE 85.5 7.17e+02 2.16e+02
+## DFOP 83.6 1.25e+11 7.64e+10
+## 
+## Representative half-life:
+## [1] 215.87
+
+
+

Example on page 6 +

+
+p6 <- nafta(NAFTA_SOP_Attachment[["p6"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p6)
+

+
+print(p6)
+
## Sums of squares:
+##       SFO      IORE      DFOP 
+## 188.45361  51.00699  42.46931 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 58.39888
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0  94.7759 7.29e-24 92.3478 97.2039
+## k_parent   0.0179 8.02e-16  0.0166  0.0194
+## sigma      3.0696 3.81e-06  2.0456  4.0936
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower    Upper
+## parent_0       97.12446 2.63e-26 95.62461 98.62431
+## k__iore_parent  0.00252 1.95e-03  0.00134  0.00472
+## N_parent        1.49587 4.07e-13  1.33896  1.65279
+## sigma           1.59698 5.05e-06  1.06169  2.13227
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0 9.66e+01 1.57e-25 95.3476 97.8979
+## k1       2.55e-02 7.33e-06  0.0233  0.0278
+## k2       3.84e-11 5.00e-01  0.0000     Inf
+## g        8.61e-01 7.55e-06  0.8314  0.8867
+## sigma    1.46e+00 6.93e-06  0.9661  1.9483
+## 
+## 
+## DTx values:
+##      DT50     DT90 DT50_rep
+## SFO  38.6 1.28e+02 3.86e+01
+## IORE 34.0 1.77e+02 5.32e+01
+## DFOP 34.1 8.50e+09 1.80e+10
+## 
+## Representative half-life:
+## [1] 53.17
+
+
+

Example on page 7 +

+
+p7 <- nafta(NAFTA_SOP_Attachment[["p7"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p7)
+

+
+print(p7)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 3661.661 3195.030 3174.145 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 3334.194
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 96.41796 4.80e-53 93.32245 99.51347
+## k_parent  0.00735 7.64e-21  0.00641  0.00843
+## sigma     7.94557 1.83e-15  6.46713  9.42401
+## 
+## $IORE
+##                Estimate Pr(>t)    Lower    Upper
+## parent_0       9.92e+01     NA 9.55e+01 1.03e+02
+## k__iore_parent 1.60e-05     NA 1.45e-07 1.77e-03
+## N_parent       2.45e+00     NA 1.35e+00 3.54e+00
+## sigma          7.42e+00     NA 6.04e+00 8.80e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower    Upper
+## parent_0 9.89e+01 9.44e-49 95.4640 102.2573
+## k1       1.81e-02 1.75e-01  0.0116   0.0281
+## k2       3.62e-10 5.00e-01  0.0000      Inf
+## g        6.06e-01 2.19e-01  0.4826   0.7178
+## sigma    7.40e+00 2.97e-15  6.0201   8.7754
+## 
+## 
+## DTx values:
+##      DT50     DT90 DT50_rep
+## SFO  94.3 3.13e+02 9.43e+01
+## IORE 96.7 1.51e+03 4.55e+02
+## DFOP 96.4 3.79e+09 1.92e+09
+## 
+## Representative half-life:
+## [1] 454.55
+
+
+
+

Examples where the representative half-life deviates from the +observed DT50 +

+
+

Example on page 8 +

+

For this dataset, the IORE fit does not converge when the default +starting values used by mkin for the IORE model are used. Therefore, a +lower value for the rate constant is used here.

+
+p8 <- nafta(NAFTA_SOP_Attachment[["p8"]], parms.ini = c(k__iore_parent = 1e-3))
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p8)
+

+
+print(p8)
+
## Sums of squares:
+##       SFO      IORE      DFOP 
+## 1996.9408  444.9237  547.5616 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 477.4924
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 88.16549 6.53e-29 83.37344 92.95754
+## k_parent  0.00803 1.67e-13  0.00674  0.00957
+## sigma     7.44786 4.17e-10  5.66209  9.23363
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower    Upper
+## parent_0       9.77e+01 7.03e-35 9.44e+01 1.01e+02
+## k__iore_parent 6.14e-05 3.20e-02 2.12e-05 1.78e-04
+## N_parent       2.27e+00 4.23e-18 2.00e+00 2.54e+00
+## sigma          3.52e+00 5.36e-10 2.67e+00 4.36e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 95.70619 8.99e-32 91.87941 99.53298
+## k1        0.02500 5.25e-04  0.01422  0.04394
+## k2        0.00273 6.84e-03  0.00125  0.00597
+## g         0.58835 2.84e-06  0.36595  0.77970
+## sigma     3.90001 6.94e-10  2.96260  4.83741
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  86.3  287     86.3
+## IORE 53.4  668    201.0
+## DFOP 55.6  517    253.0
+## 
+## Representative half-life:
+## [1] 201.03
+
+
+
+

Examples where SFO was not selected for an abiotic study +

+
+

Example on page 9, upper panel +

+
+p9a <- nafta(NAFTA_SOP_Attachment[["p9a"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p9a)
+

+
+print(p9a)
+
## Sums of squares:
+##       SFO      IORE      DFOP 
+## 839.35238  88.57064   9.93363 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 105.5678
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0  88.1933 3.06e-12 79.9447 96.4419
+## k_parent   0.0409 2.07e-07  0.0324  0.0516
+## sigma      7.2429 3.92e-05  4.4768 10.0090
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower    Upper
+## parent_0       9.89e+01 1.12e-16 9.54e+01 1.02e+02
+## k__iore_parent 1.93e-05 1.13e-01 3.49e-06 1.06e-04
+## N_parent       2.91e+00 1.45e-09 2.50e+00 3.32e+00
+## sigma          2.35e+00 5.31e-05 1.45e+00 3.26e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)  Lower  Upper
+## parent_0 9.85e+01 2.54e-20 97.390 99.672
+## k1       1.38e-01 3.52e-05  0.131  0.146
+## k2       9.02e-13 5.00e-01  0.000    Inf
+## g        6.52e-01 8.13e-06  0.642  0.661
+## sigma    7.88e-01 6.13e-02  0.481  1.095
+## 
+## 
+## DTx values:
+##      DT50     DT90 DT50_rep
+## SFO  16.9 5.63e+01 1.69e+01
+## IORE 11.6 3.37e+02 1.01e+02
+## DFOP 10.5 1.38e+12 7.68e+11
+## 
+## Representative half-life:
+## [1] 101.43
+

In this example, the residuals of the SFO indicate a lack of fit of +this model, so even if it was an abiotic experiment, the data do not +suggest a simple exponential decline.

+
+
+

Example on page 9, lower panel +

+
+p9b <- nafta(NAFTA_SOP_Attachment[["p9b"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p9b)
+

+
+print(p9b)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 35.64867 23.22334 35.64867 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 28.54188
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)  Lower   Upper
+## parent_0  94.7123 2.15e-19 93.178 96.2464
+## k_parent   0.0389 4.47e-14  0.037  0.0408
+## sigma      1.5957 1.28e-04  0.932  2.2595
+## 
+## $IORE
+##                Estimate   Pr(>t)   Lower  Upper
+## parent_0         93.863 2.32e-18 92.4565 95.269
+## k__iore_parent    0.127 1.85e-02  0.0504  0.321
+## N_parent          0.711 1.88e-05  0.4843  0.937
+## sigma             1.288 1.76e-04  0.7456  1.830
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0  94.7123 1.61e-16 93.1355 96.2891
+## k1         0.0389 1.08e-04  0.0266  0.0569
+## k2         0.0389 2.24e-04  0.0255  0.0592
+## g          0.5256 5.00e-01  0.0000  1.0000
+## sigma      1.5957 2.50e-04  0.9135  2.2779
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  17.8 59.2     17.8
+## IORE 18.4 49.2     14.8
+## DFOP 17.8 59.2     17.8
+## 
+## Representative half-life:
+## [1] 14.8
+

Here, mkin gives a longer slow DT50 for the DFOP model (17.8 days) +than PestDF (13.5 days). Presumably, this is related to the fact that +PestDF gives a negative value for the proportion of the fast degradation +which should be between 0 and 1, inclusive. This parameter is called f +in PestDF and g in mkin. In mkin, it is restricted to the interval from +0 to 1.

+
+
+

Example on page 10 +

+
+p10 <- nafta(NAFTA_SOP_Attachment[["p10"]])
+
## Warning in sqrt(diag(covar_notrans)): NaNs produced
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p10)
+

+
+print(p10)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 899.4089 336.4348 899.4089 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 413.4841
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)   Lower    Upper
+## parent_0 101.7315 6.42e-11 91.9259 111.5371
+## k_parent   0.0495 1.70e-07  0.0404   0.0607
+## sigma      8.0152 1.28e-04  4.6813  11.3491
+## 
+## $IORE
+##                Estimate   Pr(>t)  Lower   Upper
+## parent_0          96.86 3.32e-12 90.848 102.863
+## k__iore_parent     2.96 7.91e-02  0.687  12.761
+## N_parent           0.00 5.00e-01 -0.372   0.372
+## sigma              4.90 1.77e-04  2.837   6.968
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0 101.7315 1.41e-09 91.6534 111.810
+## k1         0.0495 3.04e-03  0.0188   0.131
+## k2         0.0495 4.92e-04  0.0197   0.124
+## g          0.4487      NaN  0.0000   1.000
+## sigma      8.0152 2.50e-04  4.5886  11.442
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  14.0 46.5    14.00
+## IORE 16.4 29.4     8.86
+## DFOP 14.0 46.5    14.00
+## 
+## Representative half-life:
+## [1] 8.86
+

Here, a value below N is given for the IORE model, because the data +suggests a faster decline towards the end of the experiment, which +appears physically rather unlikely in the case of a photolysis study. It +seems PestDF does not constrain N to values above zero, thus the slight +difference in IORE model parameters between PestDF and mkin.

+
+
+
+

The DT50 was not observed during the study +

+
+

Example on page 11 +

+
+p11 <- nafta(NAFTA_SOP_Attachment[["p11"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p11)
+

+
+print(p11)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 579.6805 204.7932 144.7783 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 251.6944
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 96.15820 4.83e-13 90.24934 1.02e+02
+## k_parent  0.00321 4.71e-05  0.00222 4.64e-03
+## sigma     6.43473 1.28e-04  3.75822 9.11e+00
+## 
+## $IORE
+##                Estimate Pr(>t)    Lower    Upper
+## parent_0       1.05e+02     NA 9.90e+01 1.10e+02
+## k__iore_parent 3.11e-17     NA 1.35e-20 7.18e-14
+## N_parent       8.36e+00     NA 6.62e+00 1.01e+01
+## sigma          3.82e+00     NA 2.21e+00 5.44e+00
+## 
+## $DFOP
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 1.05e+02 9.47e-13  99.9990 109.1224
+## k1       4.41e-02 5.95e-03   0.0296   0.0658
+## k2       9.94e-13 5.00e-01   0.0000      Inf
+## g        3.22e-01 1.45e-03   0.2814   0.3650
+## sigma    3.22e+00 3.52e-04   1.8410   4.5906
+## 
+## 
+## DTx values:
+##          DT50     DT90 DT50_rep
+## SFO  2.16e+02 7.18e+02 2.16e+02
+## IORE 9.73e+02 1.37e+08 4.11e+07
+## DFOP 3.07e+11 1.93e+12 6.98e+11
+## 
+## Representative half-life:
+## [1] 41148169
+

In this case, the DFOP fit reported for PestDF resulted in a negative +value for the slower rate constant, which is not possible in mkin. The +other results are in agreement.

+
+
+
+

N is less than 1 and the DFOP rate constants are like the SFO rate +constant +

+

In the following three examples, the same results are obtained with +mkin as reported for PestDF. As in the case on page 10, the N values +below 1 are deemed unrealistic and appear to be the result of an +overparameterisation.

+
+

Example on page 12, upper panel +

+
+p12a <- nafta(NAFTA_SOP_Attachment[["p12a"]])
+
## Warning in summary.mkinfit(x): Could not calculate correlation; no covariance
+## matrix
+
## Warning in sqrt(diag(covar_notrans)): NaNs produced
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p12a)
+

+
+print(p12a)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 695.4440 220.0685 695.4440 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 270.4679
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)  Lower   Upper
+## parent_0  100.521 8.75e-12 92.461 108.581
+## k_parent    0.124 3.61e-08  0.104   0.148
+## sigma       7.048 1.28e-04  4.116   9.980
+## 
+## $IORE
+##                Estimate Pr(>t) Lower Upper
+## parent_0         96.823     NA    NA    NA
+## k__iore_parent    2.436     NA    NA    NA
+## N_parent          0.263     NA    NA    NA
+## sigma             3.965     NA    NA    NA
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower   Upper
+## parent_0  100.521 2.74e-10 92.2366 108.805
+## k1          0.124 2.53e-05  0.0908   0.170
+## k2          0.124 2.52e-02  0.0456   0.339
+## g           0.793      NaN  0.0000   1.000
+## sigma       7.048 2.50e-04  4.0349  10.061
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  5.58 18.5     5.58
+## IORE 6.49 13.2     3.99
+## DFOP 5.58 18.5     5.58
+## 
+## Representative half-life:
+## [1] 3.99
+
+
+

Example on page 12, lower panel +

+
+p12b <- nafta(NAFTA_SOP_Attachment[["p12b"]])
+
## Warning in sqrt(diag(covar)): NaNs produced
+
## Warning in qt(alpha/2, rdf): NaNs produced
+
## Warning in qt(1 - alpha/2, rdf): NaNs produced
+
## Warning in pt(abs(tval), rdf, lower.tail = FALSE): NaNs produced
+
## Warning in cov2cor(ans$covar): diag(V) had non-positive or NA entries; the
+## non-finite result may be dubious
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p12b)
+

+
+print(p12b)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 58.90242 19.06353 58.90242 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 51.51756
+## 
+## Parameters:
+## $SFO
+##          Estimate  Pr(>t)   Lower    Upper
+## parent_0  97.6840 0.00039 85.9388 109.4292
+## k_parent   0.0589 0.00261  0.0431   0.0805
+## sigma      3.4323 0.04356 -1.2377   8.1023
+## 
+## $IORE
+##                Estimate Pr(>t)     Lower  Upper
+## parent_0         95.523 0.0055 74.539157 116.51
+## k__iore_parent    0.333 0.1433  0.000717 154.57
+## N_parent          0.568 0.0677 -0.989464   2.13
+## sigma             1.953 0.0975 -5.893100   9.80
+## 
+## $DFOP
+##          Estimate Pr(>t) Lower Upper
+## parent_0  97.6840    NaN   NaN   NaN
+## k1         0.0589    NaN    NA    NA
+## k2         0.0589    NaN    NA    NA
+## g          0.6473    NaN    NA    NA
+## sigma      3.4323    NaN   NaN   NaN
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  11.8 39.1    11.80
+## IORE 12.9 31.4     9.46
+## DFOP 11.8 39.1    11.80
+## 
+## Representative half-life:
+## [1] 9.46
+
+
+

Example on page 13 +

+
+p13 <- nafta(NAFTA_SOP_Attachment[["p13"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p13)
+

+
+print(p13)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 174.5971 142.3951 174.5971 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 172.131
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 92.73500 5.99e-17 89.61936 95.85065
+## k_parent  0.00258 2.42e-09  0.00223  0.00299
+## sigma     3.41172 7.07e-05  2.05455  4.76888
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower  Upper
+## parent_0        91.6016 6.34e-16 88.53086 94.672
+## k__iore_parent   0.0396 2.36e-01  0.00207  0.759
+## N_parent         0.3541 1.46e-01 -0.35153  1.060
+## sigma            3.0811 9.64e-05  1.84296  4.319
+## 
+## $DFOP
+##          Estimate Pr(>t)    Lower    Upper
+## parent_0 92.73500     NA 8.95e+01 95.92118
+## k1        0.00258     NA 4.18e-04  0.01592
+## k2        0.00258     NA 1.75e-03  0.00381
+## g         0.16452     NA 0.00e+00  1.00000
+## sigma     3.41172     NA 2.02e+00  4.79960
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO   269  892      269
+## IORE  261  560      169
+## DFOP  269  892      269
+## 
+## Representative half-life:
+## [1] 168.51
+
+
+
+

DT50 not observed in the study and DFOP problems in PestDF +

+
+p14 <- nafta(NAFTA_SOP_Attachment[["p14"]])
+
## Warning in sqrt(diag(covar)): NaNs produced
+
## Warning in cov2cor(ans$covar): diag(V) had non-positive or NA entries; the
+## non-finite result may be dubious
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p14)
+

+
+print(p14)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 48.43249 28.67746 27.26248 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 32.83337
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 99.47124 2.06e-30 98.42254 1.01e+02
+## k_parent  0.00279 3.75e-15  0.00256 3.04e-03
+## sigma     1.55616 3.81e-06  1.03704 2.08e+00
+## 
+## $IORE
+##                Estimate Pr(>t) Lower Upper
+## parent_0       1.00e+02     NA   NaN   NaN
+## k__iore_parent 9.44e-08     NA   NaN   NaN
+## N_parent       3.31e+00     NA   NaN   NaN
+## sigma          1.20e+00     NA 0.796   1.6
+## 
+## $DFOP
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 1.00e+02 2.96e-28 99.40280 101.2768
+## k1       9.53e-03 1.20e-01  0.00638   0.0143
+## k2       5.21e-12 5.00e-01  0.00000      Inf
+## g        3.98e-01 2.19e-01  0.30481   0.4998
+## sigma    1.17e+00 7.68e-06  0.77406   1.5610
+## 
+## 
+## DTx values:
+##          DT50     DT90 DT50_rep
+## SFO  2.48e+02 8.25e+02 2.48e+02
+## IORE 4.34e+02 2.22e+04 6.70e+03
+## DFOP 3.55e+10 3.44e+11 1.33e+11
+## 
+## Representative half-life:
+## [1] 6697.44
+

The slower rate constant reported by PestDF is negative, which is not +physically realistic, and not possible in mkin. The other fits give the +same results in mkin and PestDF.

+
+
+

N is less than 1 and DFOP fraction parameter is below zero +

+
+p15a <- nafta(NAFTA_SOP_Attachment[["p15a"]])
+
## Warning in sqrt(diag(covar)): NaNs produced
+
## Warning in cov2cor(ans$covar): diag(V) had non-positive or NA entries; the
+## non-finite result may be dubious
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p15a)
+

+
+print(p15a)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 245.5248 135.0132 245.5248 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 165.9335
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower   Upper
+## parent_0 97.96751 2.00e-15 94.32049 101.615
+## k_parent  0.00952 4.93e-09  0.00824   0.011
+## sigma     4.18778 1.28e-04  2.44588   5.930
+## 
+## $IORE
+##                Estimate   Pr(>t)  Lower  Upper
+## parent_0         95.874 2.94e-15 92.937 98.811
+## k__iore_parent    0.629 2.11e-01  0.044  8.982
+## N_parent          0.000 5.00e-01 -0.642  0.642
+## sigma             3.105 1.78e-04  1.795  4.416
+## 
+## $DFOP
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 97.96751 2.85e-13 94.21913 101.7159
+## k1        0.00952 6.28e-02  0.00260   0.0349
+## k2        0.00952 1.27e-04  0.00652   0.0139
+## g         0.21241 5.00e-01       NA       NA
+## sigma     4.18778 2.50e-04  2.39747   5.9781
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  72.8  242     72.8
+## IORE 76.3  137     41.3
+## DFOP 72.8  242     72.8
+## 
+## Representative half-life:
+## [1] 41.33
+
+p15b <- nafta(NAFTA_SOP_Attachment[["p15b"]])
+
## Warning in summary.mkinfit(x): Could not calculate correlation; no covariance
+## matrix
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The half-life obtained from the IORE model may be used
+
+plot(p15b)
+

+
+print(p15b)
+
## Sums of squares:
+##       SFO      IORE      DFOP 
+## 106.91629  68.55574 106.91629 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 84.25618
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)    Lower    Upper
+## parent_0 1.01e+02 3.06e-17 98.31594 1.03e+02
+## k_parent 4.86e-03 2.48e-10  0.00435 5.42e-03
+## sigma    2.76e+00 1.28e-04  1.61402 3.91e+00
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower  Upper
+## parent_0          99.83 1.81e-16 97.51348 102.14
+## k__iore_parent     0.38 3.22e-01  0.00352  41.05
+## N_parent           0.00 5.00e-01 -1.07696   1.08
+## sigma              2.21 2.57e-04  1.23245   3.19
+## 
+## $DFOP
+##          Estimate Pr(>t) Lower Upper
+## parent_0 1.01e+02     NA    NA    NA
+## k1       4.86e-03     NA    NA    NA
+## k2       4.86e-03     NA    NA    NA
+## g        1.88e-01     NA    NA    NA
+## sigma    2.76e+00     NA    NA    NA
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO   143  474    143.0
+## IORE  131  236     71.2
+## DFOP  143  474    143.0
+## 
+## Representative half-life:
+## [1] 71.18
+

In mkin, only the IORE fit is affected (deemed unrealistic), as the +fraction parameter of the DFOP model is restricted to the interval +between 0 and 1 in mkin. The SFO fits give the same results for both +mkin and PestDF.

+
+
+

The DFOP fraction parameter is greater than 1 +

+
+p16 <- nafta(NAFTA_SOP_Attachment[["p16"]])
+
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
+
## The representative half-life of the IORE model is longer than the one corresponding
+
## to the terminal degradation rate found with the DFOP model.
+
## The representative half-life obtained from the DFOP model may be used
+
+plot(p16)
+

+
+print(p16)
+
## Sums of squares:
+##      SFO     IORE     DFOP 
+## 3831.804 2062.008 1550.980 
+## 
+## Critical sum of squares for checking the SFO model:
+## [1] 2247.348
+## 
+## Parameters:
+## $SFO
+##          Estimate   Pr(>t)  Lower Upper
+## parent_0   71.953 2.33e-13 60.509 83.40
+## k_parent    0.159 4.86e-05  0.102  0.25
+## sigma      11.302 1.25e-08  8.308 14.30
+## 
+## $IORE
+##                Estimate   Pr(>t)    Lower    Upper
+## parent_0       8.74e+01 2.48e-16 7.72e+01 97.52972
+## k__iore_parent 4.55e-04 2.16e-01 3.48e-05  0.00595
+## N_parent       2.70e+00 1.21e-08 1.99e+00  3.40046
+## sigma          8.29e+00 1.61e-08 6.09e+00 10.49062
+## 
+## $DFOP
+##          Estimate   Pr(>t)   Lower  Upper
+## parent_0  88.5333 7.40e-18 79.9836 97.083
+## k1        18.8461 5.00e-01  0.0000    Inf
+## k2         0.0776 1.41e-05  0.0518  0.116
+## g          0.4733 1.41e-09  0.3674  0.582
+## sigma      7.1902 2.11e-08  5.2785  9.102
+## 
+## 
+## DTx values:
+##      DT50 DT90 DT50_rep
+## SFO  4.35 14.4     4.35
+## IORE 1.48 32.1     9.67
+## DFOP 0.67 21.4     8.93
+## 
+## Representative half-life:
+## [1] 8.93
+

In PestDF, the DFOP fit seems to have stuck in a local minimum, as +mkin finds a solution with a much lower +χ2\chi^2 +error level. As the half-life from the slower rate constant of the DFOP +model is larger than the IORE derived half-life, the NAFTA +recommendation obtained with mkin is to use the DFOP representative +half-life of 8.9 days.

+
+
+

Conclusions +

+

The results obtained with mkin deviate from the results obtained with +PestDF either in cases where one of the interpretive rules would apply, +i.e. the IORE parameter N is less than one or the DFOP k values obtained +with PestDF are equal to the SFO k values, or in cases where the DFOP +model did not converge, which often lead to negative rate constants +returned by PestDF.

+

Therefore, mkin appears to suitable for kinetic evaluations according +to the NAFTA guidance.

+
+
+

References +

+
+
+US EPA. 2015. “Standard Operating Procedure for Using the NAFTA +Guidance to Calculate Representative Half-Life Values and Characterizing +Pesticide Degradation.” https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance. +
+
+
+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p10-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p10-1.png new file mode 100644 index 00000000..1d4a25e0 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p10-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p11-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p11-1.png new file mode 100644 index 00000000..71fc4699 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p11-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12a-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12a-1.png new file mode 100644 index 00000000..a1d3a084 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12a-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12b-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12b-1.png new file mode 100644 index 00000000..1a6fdd03 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p12b-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p13-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p13-1.png new file mode 100644 index 00000000..f9b9f637 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p13-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png new file mode 100644 index 00000000..9f7b0cc5 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15a-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15a-1.png new file mode 100644 index 00000000..aa55169e Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15a-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15b-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15b-1.png new file mode 100644 index 00000000..d17c7aae Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p15b-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p16-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p16-1.png new file mode 100644 index 00000000..75ac7e5b Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p16-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5a-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5a-1.png new file mode 100644 index 00000000..12a62954 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5a-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5b-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5b-1.png new file mode 100644 index 00000000..9e38e696 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p5b-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p6-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p6-1.png new file mode 100644 index 00000000..e6e3abbe Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p6-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p7-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p7-1.png new file mode 100644 index 00000000..7c5d4bab Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p7-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p8-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p8-1.png new file mode 100644 index 00000000..a1e3bf25 Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p8-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9a-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9a-1.png new file mode 100644 index 00000000..c247fd4e Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9a-1.png differ diff --git a/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png new file mode 100644 index 00000000..99d593fc Binary files /dev/null and b/docs/dev/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png differ diff --git a/docs/dev/articles/web_only/benchmarks.html b/docs/dev/articles/web_only/benchmarks.html new file mode 100644 index 00000000..bd6a443a --- /dev/null +++ b/docs/dev/articles/web_only/benchmarks.html @@ -0,0 +1,1133 @@ + + + + + + + +Benchmark timings for mkin • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Benchmark timings for mkin

+

Johannes +Ranke

+ +

Last change 17 February 2023 +(rebuilt 2025-02-14)

+ + Source: vignettes/web_only/benchmarks.rmd +
benchmarks.rmd
+
+ + + +

Each system is characterized by the operating system type, the CPU +type, the mkin version, and, as in June 2022 the current R version lead +to worse performance, the R version. A compiler was available, so if no +analytical solution was available, compiled ODE models are used.

+

Every fit is only performed once, so the accuracy of the benchmarks +is limited.

+

The following wrapper function for mmkin is used because +the way the error model is specified was changed in mkin version +0.9.49.1.

+
+if (packageVersion("mkin") > "0.9.48.1") {
+  mmkin_bench <- function(models, datasets, error_model = "const") {
+    mmkin(models, datasets, error_model = error_model, cores = 1, quiet = TRUE)
+  }
+} else {
+  mmkin_bench <- function(models, datasets, error_model = NULL) {
+    mmkin(models, datasets, reweight.method = error_model, cores = 1, quiet = TRUE)
+  }
+}
+
+

Test cases +

+

Parent only:

+
+FOCUS_C <- FOCUS_2006_C
+FOCUS_D <- subset(FOCUS_2006_D, value != 0)
+parent_datasets <- list(FOCUS_C, FOCUS_D)
+
+
+t1 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), parent_datasets))[["elapsed"]]
+t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), parent_datasets,
+    error_model = "tc"))[["elapsed"]]
+

One metabolite:

+
+SFO_SFO <- mkinmod(
+  parent = mkinsub("SFO", "m1"),
+  m1 = mkinsub("SFO"))
+FOMC_SFO <- mkinmod(
+  parent = mkinsub("FOMC", "m1"),
+  m1 = mkinsub("SFO"))
+DFOP_SFO <- mkinmod(
+  parent = mkinsub("FOMC", "m1"), # erroneously used FOMC twice, not fixed for consistency
+  m1 = mkinsub("SFO"))
+t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D)))[["elapsed"]]
+t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D),
+    error_model = "tc"))[["elapsed"]]
+t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D),
+    error_model = "obs"))[["elapsed"]]
+

Two metabolites, synthetic data:

+
+m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"),
+                           M1 = mkinsub("SFO", "M2"),
+                           M2 = mkinsub("SFO"),
+                           use_of_ff = "max", quiet = TRUE)
+
+m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
+                           M1 = mkinsub("SFO"),
+                           M2 = mkinsub("SFO"),
+                           use_of_ff = "max", quiet = TRUE)
+
+SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
+
+DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
+
+t6 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a)))[["elapsed"]]
+t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))[["elapsed"]]
+
+t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
+    error_model = "tc"))[["elapsed"]]
+t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
+    error_model = "tc"))[["elapsed"]]
+
+t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a),
+    error_model = "obs"))[["elapsed"]]
+t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c),
+    error_model = "obs"))[["elapsed"]]
+
+
+

Results +

+

Benchmarks for all available error models are shown. They are +intended for improving mkin, not for comparing CPUs or operating +systems. All trademarks belong to their respective owners.

+
+

Parent only +

+

Constant variance (t1) and two-component error model (t2) for four +models fitted to two datasets, i.e. eight fits for each test.

+ ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
OSCPURmkint1t2
LinuxRyzen 7 1700NA0.9.48.13.61011.019
LinuxRyzen 7 1700NA0.9.49.18.18422.889
LinuxRyzen 7 1700NA0.9.49.27.06412.558
LinuxRyzen 7 1700NA0.9.49.37.29621.239
LinuxRyzen 7 1700NA0.9.49.45.93620.545
LinuxRyzen 7 1700NA0.9.50.21.7143.971
LinuxRyzen 7 1700NA0.9.50.31.7524.156
LinuxRyzen 7 1700NA0.9.50.41.7863.729
LinuxRyzen 7 1700NA1.0.31.8813.504
LinuxRyzen 7 1700NA1.0.41.8673.450
LinuxRyzen 7 17004.1.31.1.01.7913.289
LinuxRyzen 7 17004.2.11.1.01.8423.453
Linuxi7-4710MQ4.2.11.1.01.9594.116
Linuxi7-4710MQ4.1.31.1.01.8773.906
Linuxi7-4710MQ4.2.11.1.11.6443.172
LinuxRyzen 7 17004.2.11.1.11.7703.377
LinuxRyzen 7 17004.2.11.1.21.9573.633
LinuxRyzen 7 17004.2.21.2.02.1403.774
LinuxRyzen 7 17004.2.21.2.22.1873.851
LinuxRyzen 9 7950X4.2.21.2.01.2881.794
LinuxRyzen 9 7950X4.2.21.2.21.2761.804
LinuxRyzen 9 7950X4.2.21.2.31.3701.883
LinuxRyzen 9 7950X4.2.31.2.31.4061.948
LinuxRyzen 9 7950X4.3.01.2.41.3861.960
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.52.3693.632
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.62.8564.960
LinuxRyzen 9 7950X4.3.21.2.61.4082.041
LinuxRyzen 9 7950X4.4.21.2.91.3231.925
LinuxRyzen 9 7950X4.4.21.2.101.3711.980
+
+
+

One metabolite +

+

Constant variance (t3), two-component error model (t4), and variance +by variable (t5) for three models fitted to one dataset, i.e. three fits +for each test.

+ +++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
OSCPURmkint3t4t5
LinuxRyzen 7 1700NA0.9.48.13.76414.3479.495
LinuxRyzen 7 1700NA0.9.49.14.64913.7896.395
LinuxRyzen 7 1700NA0.9.49.24.7868.4615.675
LinuxRyzen 7 1700NA0.9.49.34.51013.8057.386
LinuxRyzen 7 1700NA0.9.49.44.44615.3356.002
LinuxRyzen 7 1700NA0.9.50.21.4026.1742.764
LinuxRyzen 7 1700NA0.9.50.31.4306.6152.878
LinuxRyzen 7 1700NA0.9.50.41.3977.2512.810
LinuxRyzen 7 1700NA1.0.31.4306.3442.798
LinuxRyzen 7 1700NA1.0.41.4156.3642.820
LinuxRyzen 7 17004.1.31.1.01.3106.2792.681
LinuxRyzen 7 17004.2.11.1.03.80221.2478.461
Linuxi7-4710MQ4.2.11.1.03.33419.5217.565
Linuxi7-4710MQ4.1.31.1.01.5788.0583.339
Linuxi7-4710MQ4.2.11.1.11.2305.8392.444
LinuxRyzen 7 17004.2.11.1.11.3085.7582.558
LinuxRyzen 7 17004.2.11.1.21.5036.1472.803
LinuxRyzen 7 17004.2.21.2.01.5546.1932.843
LinuxRyzen 7 17004.2.21.2.21.5856.3353.003
LinuxRyzen 9 7950X4.2.21.2.00.7922.3781.245
LinuxRyzen 9 7950X4.2.21.2.20.7842.3551.233
LinuxRyzen 9 7950X4.2.21.2.30.7702.0111.123
LinuxRyzen 9 7950X4.2.31.2.30.7932.1091.178
LinuxRyzen 9 7950X4.3.01.2.40.7792.0801.106
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.51.8235.5552.404
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.61.7615.4052.462
LinuxRyzen 9 7950X4.3.21.2.60.7952.2281.178
LinuxRyzen 9 7950X4.4.21.2.90.7542.1531.139
LinuxRyzen 9 7950X4.4.21.2.100.7622.1791.131
+
+
+

Two metabolites +

+

Constant variance (t6 and t7), two-component error model (t8 and t9), +and variance by variable (t10 and t11) for one model fitted to one +dataset, i.e. one fit for each test.

+ ++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
OSCPURmkint6t7t8t9t10t11
LinuxRyzen 7 1700NA0.9.48.12.6234.5877.52516.6218.57631.267
LinuxRyzen 7 1700NA0.9.49.12.5424.1284.6328.1713.6765.636
LinuxRyzen 7 1700NA0.9.49.22.7234.4784.8627.6183.5795.574
LinuxRyzen 7 1700NA0.9.49.32.6434.3747.02011.1245.3887.365
LinuxRyzen 7 1700NA0.9.49.42.6354.2594.7377.7633.4275.626
LinuxRyzen 7 1700NA0.9.50.20.7771.2361.3322.8722.0692.987
LinuxRyzen 7 1700NA0.9.50.30.8581.2641.3332.9842.1133.073
LinuxRyzen 7 1700NA0.9.50.40.7831.2821.4863.8151.9583.105
LinuxRyzen 7 1700NA1.0.30.7631.2441.4573.0541.9232.839
LinuxRyzen 7 1700NA1.0.40.7851.2521.4663.0911.9362.826
LinuxRyzen 7 17004.1.31.1.00.7441.2271.2883.5531.8952.738
LinuxRyzen 7 17004.2.11.1.03.0184.1655.03610.8446.6239.722
Linuxi7-4710MQ4.2.11.1.02.5223.7924.14311.2685.9358.728
Linuxi7-4710MQ4.1.31.1.00.9071.5351.5894.5442.3023.463
Linuxi7-4710MQ4.2.11.1.10.6781.0951.1493.2471.6582.472
LinuxRyzen 7 17004.2.11.1.10.6961.1241.3212.7861.7442.566
LinuxRyzen 7 17004.2.11.1.20.8611.2951.5073.1021.9612.852
LinuxRyzen 7 17004.2.21.2.00.9131.3451.5393.0111.9872.802
LinuxRyzen 7 17004.2.21.2.20.9351.3811.5513.2091.9763.013
LinuxRyzen 9 7950X4.2.21.2.00.4450.5910.6601.1900.8141.100
LinuxRyzen 9 7950X4.2.21.2.20.4430.5860.6611.1760.8031.097
LinuxRyzen 9 7950X4.2.21.2.30.4180.5300.5911.0060.7160.949
LinuxRyzen 9 7950X4.2.31.2.30.4320.5490.6091.0520.7430.989
LinuxRyzen 9 7950X4.3.01.2.40.4100.5260.5531.2490.7120.948
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.50.7981.0961.2173.1731.6342.271
LinuxIntel(R) Xeon(R) Gold 6134 CPU @ 3.20GHz4.3.11.2.60.8131.1361.2203.1141.5982.255
LinuxRyzen 9 7950X4.3.21.2.60.4390.5570.5851.3380.7490.999
LinuxRyzen 9 7950X4.4.21.2.90.4240.5340.5601.2980.7350.981
LinuxRyzen 9 7950X4.4.21.2.100.4260.5330.5651.2960.7280.978
+
+
+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/compiled_models.html b/docs/dev/articles/web_only/compiled_models.html new file mode 100644 index 00000000..113f7b67 --- /dev/null +++ b/docs/dev/articles/web_only/compiled_models.html @@ -0,0 +1,231 @@ + + + + + + + +Performance benefit by using compiled model definitions in mkin • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Performance benefit by using compiled model definitions in mkin

+

Johannes +Ranke

+ +

2025-02-14

+ + Source: vignettes/web_only/compiled_models.rmd +
compiled_models.rmd
+
+ + + +
+

How to benefit from compiled models +

+

When using an mkin version equal to or greater than 0.9-36 and a C +compiler is available, you will see a message that the model is being +compiled from autogenerated C code when defining a model using mkinmod. +Starting from version 0.9.49.9, the mkinmod() function +checks for presence of a compiler using

+
+pkgbuild::has_compiler()
+

In previous versions, it used Sys.which("gcc") for this +check.

+

On Linux, you need to have the essential build tools like make and +gcc or clang installed. On Debian based linux distributions, these will +be pulled in by installing the build-essential package.

+

On MacOS, which I do not use personally, I have had reports that a +compiler is available by default.

+

On Windows, you need to install Rtools and have the path to its bin +directory in your PATH variable. You do not need to modify the PATH +variable when installing Rtools. Instead, I would recommend to put the +line

+
+Sys.setenv(PATH = paste("C:/Rtools/bin", Sys.getenv("PATH"), sep=";"))
+

into your .Rprofile startup file. This is just a text file with some +R code that is executed when your R session starts. It has to be named +.Rprofile and has to be located in your home directory, which will +generally be your Documents folder. You can check the location of the +home directory used by R by issuing

+
+Sys.getenv("HOME")
+
+
+

Comparison with other solution methods +

+

First, we build a simple degradation model for a parent compound with +one metabolite, and we remove zero values from the dataset.

+
+library("mkin", quietly = TRUE)
+SFO_SFO <- mkinmod(
+  parent = mkinsub("SFO", "m1"),
+  m1 = mkinsub("SFO"))
+
## Temporary DLL for differentials generated and loaded
+
+FOCUS_D <- subset(FOCUS_2006_D, value != 0)
+

We can compare the performance of the Eigenvalue based solution +against the compiled version and the R implementation of the +differential equations using the benchmark package. In the output of +below code, the warnings about zero being removed from the FOCUS D +dataset are suppressed. Since mkin version 0.9.49.11, an analytical +solution is also implemented, which is included in the tests below.

+
+if (require(rbenchmark)) {
+  b.1 <- benchmark(
+    "deSolve, not compiled" = mkinfit(SFO_SFO, FOCUS_D,
+       solution_type = "deSolve",
+       use_compiled = FALSE, quiet = TRUE),
+    "Eigenvalue based" = mkinfit(SFO_SFO, FOCUS_D,
+       solution_type = "eigen", quiet = TRUE),
+    "deSolve, compiled" = mkinfit(SFO_SFO, FOCUS_D,
+       solution_type = "deSolve", quiet = TRUE),
+    "analytical" = mkinfit(SFO_SFO, FOCUS_D,
+       solution_type = "analytical",
+       use_compiled = FALSE, quiet = TRUE),
+    replications = 1, order = "relative",
+    columns = c("test", "replications", "relative", "elapsed"))
+  print(b.1)
+} else {
+  print("R package rbenchmark is not available")
+}
+
##                    test replications relative elapsed
+## 4            analytical            1    1.000   0.102
+## 3     deSolve, compiled            1    1.324   0.135
+## 2      Eigenvalue based            1    1.706   0.174
+## 1 deSolve, not compiled            1   22.627   2.308
+

We see that using the compiled model is by more than a factor of 10 +faster than using deSolve without compiled code.

+
+
+

Model without analytical solution +

+

This evaluation is also taken from the example section of mkinfit. No +analytical solution is available for this system, and now Eigenvalue +based solution is possible, so only deSolve using with or without +compiled code is available.

+
+if (require(rbenchmark)) {
+  FOMC_SFO <- mkinmod(
+    parent = mkinsub("FOMC", "m1"),
+    m1 = mkinsub( "SFO"))
+
+  b.2 <- benchmark(
+    "deSolve, not compiled" = mkinfit(FOMC_SFO, FOCUS_D,
+                                      use_compiled = FALSE, quiet = TRUE),
+    "deSolve, compiled" = mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE),
+    replications = 1, order = "relative",
+    columns = c("test", "replications", "relative", "elapsed"))
+  print(b.2)
+  factor_FOMC_SFO <- round(b.2["1", "relative"])
+} else {
+  factor_FOMC_SFO <- NA
+  print("R package benchmark is not available")
+}
+
## Temporary DLL for differentials generated and loaded
+
##                    test replications relative elapsed
+## 2     deSolve, compiled            1    1.000   0.170
+## 1 deSolve, not compiled            1   23.865   4.057
+

Here we get a performance benefit of a factor of 24 using the version +of the differential equation model compiled from C code!

+

This vignette was built with mkin 1.2.10 on

+
## R version 4.4.2 (2024-10-31)
+## Platform: x86_64-pc-linux-gnu
+## Running under: Debian GNU/Linux 12 (bookworm)
+
## CPU model: AMD Ryzen 9 7950X 16-Core Processor
+
+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/dimethenamid_2018.html b/docs/dev/articles/web_only/dimethenamid_2018.html new file mode 100644 index 00000000..572f253d --- /dev/null +++ b/docs/dev/articles/web_only/dimethenamid_2018.html @@ -0,0 +1,651 @@ + + + + + + + +Example evaluations of the dimethenamid data from 2018 • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Example evaluations of the dimethenamid data from 2018

+

Johannes +Ranke

+ +

Last change 1 July 2022, +built on 14 Feb 2025

+ + Source: vignettes/web_only/dimethenamid_2018.rmd +
dimethenamid_2018.rmd
+
+ + + +

Wissenschaftlicher Berater, Kronacher +Str. 12, 79639 Grenzach-Wyhlen, Germany

+
+

Introduction +

+

A first analysis of the data analysed here was presented in a recent +journal article on nonlinear mixed-effects models in degradation +kinetics (Ranke et al. 2021). That +analysis was based on the nlme package and a development +version of the saemix package that was unpublished at the +time. Meanwhile, version 3.0 of the saemix package is +available from the CRAN repository. Also, it turned out that there was +an error in the handling of the Borstel data in the mkin package at the +time, leading to the duplication of a few data points from that soil. +The dataset in the mkin package has been corrected, and the interface to +saemix in the mkin package has been updated to use the +released version.

+

This vignette is intended to present an up to date analysis of the +data, using the corrected dataset and released versions of +mkin and saemix.

+
+
+

Data +

+

Residue data forming the basis for the endpoints derived in the +conclusion on the peer review of the pesticide risk assessment of +dimethenamid-P published by the European Food Safety Authority (EFSA) in +2018 (EFSA 2018) were transcribed from the +risk assessment report (Rapporteur Member State +Germany, Co-Rapporteur Member State Bulgaria 2018) which can be +downloaded from the Open EFSA repository https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716.

+

The data are available +in the mkin package. The following code (hidden by default, please +use the button to the right to show it) treats the data available for +the racemic mixture dimethenamid (DMTA) and its enantiomer +dimethenamid-P (DMTAP) in the same way, as no difference between their +degradation behaviour was identified in the EU risk assessment. The +observation times of each dataset are multiplied with the corresponding +normalisation factor also available in the dataset, in order to make it +possible to describe all datasets with a single set of parameters.

+

Also, datasets observed in the same soil are merged, resulting in +dimethenamid (DMTA) data from six soils.

+
+library(mkin, quietly = TRUE)
+dmta_ds <- lapply(1:7, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+  ds_i
+})
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL
+
+
+

Parent degradation +

+

We evaluate the observed degradation of the parent compound using +simple exponential decline (SFO) and biexponential decline (DFOP), using +constant variance (const) and a two-component variance (tc) as error +models.

+
+

Separate evaluations +

+

As a first step, to get a visual impression of the fit of the +different models, we do separate evaluations for each soil using the +mmkin function from the mkin package:

+
+f_parent_mkin_const <- mmkin(c("SFO", "DFOP"), dmta_ds,
+  error_model = "const", quiet = TRUE)
+f_parent_mkin_tc <- mmkin(c("SFO", "DFOP"), dmta_ds,
+  error_model = "tc", quiet = TRUE)
+

The plot of the individual SFO fits shown below suggests that at +least in some datasets the degradation slows down towards later time +points, and that the scatter of the residuals error is smaller for +smaller values (panel to the right):

+
+plot(mixed(f_parent_mkin_const["SFO", ]))
+

+

Using biexponential decline (DFOP) results in a slightly more random +scatter of the residuals:

+
+plot(mixed(f_parent_mkin_const["DFOP", ]))
+

+

The population curve (bold line) in the above plot results from +taking the mean of the individual transformed parameters, i.e. of log k1 +and log k2, as well as of the logit of the g parameter of the DFOP +model). Here, this procedure does not result in parameters that +represent the degradation well, because in some datasets the fitted +value for k2 is extremely close to zero, leading to a log k2 value that +dominates the average. This is alleviated if only rate constants that +pass the t-test for significant difference from zero (on the +untransformed scale) are considered in the averaging:

+
+plot(mixed(f_parent_mkin_const["DFOP", ]), test_log_parms = TRUE)
+

+

While this is visually much more satisfactory, such an average +procedure could introduce a bias, as not all results from the individual +fits enter the population curve with the same weight. This is where +nonlinear mixed-effects models can help out by treating all datasets +with equally by fitting a parameter distribution model together with the +degradation model and the error model (see below).

+

The remaining trend of the residuals to be higher for higher +predicted residues is reduced by using the two-component error +model:

+
+plot(mixed(f_parent_mkin_tc["DFOP", ]), test_log_parms = TRUE)
+

+

However, note that in the case of using this error model, the fits to +the Flaach and BBA 2.3 datasets appear to be ill-defined, indicated by +the fact that they did not converge:

+
+print(f_parent_mkin_tc["DFOP", ])
+
<mmkin> object
+Status of individual fits:
+
+      dataset
+model  Calke Borstel Flaach BBA 2.2 BBA 2.3 Elliot
+  DFOP OK    OK      OK     OK      C       OK    
+
+C: Optimisation did not converge:
+iteration limit reached without convergence (10)
+OK: No warnings
+
+
+

Nonlinear mixed-effects models +

+

Instead of taking a model selection decision for each of the +individual fits, we fit nonlinear mixed-effects models (using different +fitting algorithms as implemented in different packages) and do model +selection using all available data at the same time. In order to make +sure that these decisions are not unduly influenced by the type of +algorithm used, by implementation details or by the use of wrong control +parameters, we compare the model selection results obtained with +different R packages, with different algorithms and checking control +parameters.

+
+

nlme +

+

The nlme package was the first R extension providing facilities to +fit nonlinear mixed-effects models. We would like to do model selection +from all four combinations of degradation models and error models based +on the AIC. However, fitting the DFOP model with constant variance and +using default control parameters results in an error, signalling that +the maximum number of 50 iterations was reached, potentially indicating +overparameterisation. Nevertheless, the algorithm converges when the +two-component error model is used in combination with the DFOP model. +This can be explained by the fact that the smaller residues observed at +later sampling times get more weight when using the two-component error +model which will counteract the tendency of the algorithm to try +parameter combinations unsuitable for fitting these data.

+
+library(nlme)
+f_parent_nlme_sfo_const <- nlme(f_parent_mkin_const["SFO", ])
+# f_parent_nlme_dfop_const <- nlme(f_parent_mkin_const["DFOP", ])
+f_parent_nlme_sfo_tc <- nlme(f_parent_mkin_tc["SFO", ])
+f_parent_nlme_dfop_tc <- nlme(f_parent_mkin_tc["DFOP", ])
+

Note that a certain degree of overparameterisation is also indicated +by a warning obtained when fitting DFOP with the two-component error +model (‘false convergence’ in the ‘LME step’ in iteration 3). However, +as this warning does not occur in later iterations, and specifically not +in the last of the 5 iterations, we can ignore this warning.

+

The model comparison function of the nlme package can directly be +applied to these fits showing a much lower AIC for the DFOP model fitted +with the two-component error model. Also, the likelihood ratio test +indicates that this difference is significant as the p-value is below +0.0001.

+
+anova(
+  f_parent_nlme_sfo_const, f_parent_nlme_sfo_tc, f_parent_nlme_dfop_tc
+)
+
                        Model df    AIC    BIC  logLik   Test L.Ratio p-value
+f_parent_nlme_sfo_const     1  5 796.60 811.82 -393.30                       
+f_parent_nlme_sfo_tc        2  6 798.60 816.86 -393.30 1 vs 2    0.00   0.998
+f_parent_nlme_dfop_tc       3 10 671.91 702.34 -325.95 2 vs 3  134.69  <.0001
+

In addition to these fits, attempts were also made to include +correlations between random effects by using the log Cholesky +parameterisation of the matrix specifying them. The code used for these +attempts can be made visible below.

+
+f_parent_nlme_sfo_const_logchol <- nlme(f_parent_mkin_const["SFO", ],
+  random = nlme::pdLogChol(list(DMTA_0 ~ 1, log_k_DMTA ~ 1)))
+anova(f_parent_nlme_sfo_const, f_parent_nlme_sfo_const_logchol)
+f_parent_nlme_sfo_tc_logchol <- nlme(f_parent_mkin_tc["SFO", ],
+  random = nlme::pdLogChol(list(DMTA_0 ~ 1, log_k_DMTA ~ 1)))
+anova(f_parent_nlme_sfo_tc, f_parent_nlme_sfo_tc_logchol)
+f_parent_nlme_dfop_tc_logchol <- nlme(f_parent_mkin_const["DFOP", ],
+  random = nlme::pdLogChol(list(DMTA_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1)))
+anova(f_parent_nlme_dfop_tc, f_parent_nlme_dfop_tc_logchol)
+

While the SFO variants converge fast, the additional parameters +introduced by this lead to convergence warnings for the DFOP model. The +model comparison clearly show that adding correlations between random +effects does not improve the fits.

+

The selected model (DFOP with two-component error) fitted to the data +assuming no correlations between random effects is shown below.

+
+plot(f_parent_nlme_dfop_tc)
+

+
+
+

saemix +

+

The saemix package provided the first Open Source implementation of +the Stochastic Approximation to the Expectation Maximisation (SAEM) +algorithm. SAEM fits of degradation models can be conveniently performed +using an interface to the saemix package available in current +development versions of the mkin package.

+

The corresponding SAEM fits of the four combinations of degradation +and error models are fitted below. As there is no convergence criterion +implemented in the saemix package, the convergence plots need to be +manually checked for every fit. We define control settings that work +well for all the parent data fits shown in this vignette.

+
+library(saemix)
+saemix_control <- saemixControl(nbiter.saemix = c(800, 300), nb.chains = 15,
+    print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)
+saemix_control_moreiter <- saemixControl(nbiter.saemix = c(1600, 300), nb.chains = 15,
+    print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)
+saemix_control_10k <- saemixControl(nbiter.saemix = c(10000, 300), nb.chains = 15,
+    print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)
+

The convergence plot for the SFO model using constant variance is +shown below.

+
+f_parent_saemix_sfo_const <- mkin::saem(f_parent_mkin_const["SFO", ], quiet = TRUE,
+  control = saemix_control, transformations = "saemix")
+plot(f_parent_saemix_sfo_const$so, plot.type = "convergence")
+

+

Obviously the selected number of iterations is sufficient to reach +convergence. This can also be said for the SFO fit using the +two-component error model.

+
+f_parent_saemix_sfo_tc <- mkin::saem(f_parent_mkin_tc["SFO", ], quiet = TRUE,
+  control = saemix_control, transformations = "saemix")
+plot(f_parent_saemix_sfo_tc$so, plot.type = "convergence")
+

+

When fitting the DFOP model with constant variance (see below), +parameter convergence is not as unambiguous.

+
+f_parent_saemix_dfop_const <- mkin::saem(f_parent_mkin_const["DFOP", ], quiet = TRUE,
+  control = saemix_control, transformations = "saemix")
+plot(f_parent_saemix_dfop_const$so, plot.type = "convergence")
+

+
+print(f_parent_saemix_dfop_const)
+
Kinetic nonlinear mixed-effects model fit by SAEM
+Structural model:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+           time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+           * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Likelihood computed by importance sampling
+  AIC BIC logLik
+  706 704   -344
+
+Fitted parameters:
+          estimate    lower   upper
+DMTA_0    97.99583 96.50079 99.4909
+k1         0.06377  0.03432  0.0932
+k2         0.00848  0.00444  0.0125
+g          0.95701  0.91313  1.0009
+a.1        1.82141  1.60516  2.0377
+SD.DMTA_0  1.64787  0.45729  2.8384
+SD.k1      0.57439  0.24731  0.9015
+SD.k2      0.03296 -2.50524  2.5712
+SD.g       1.10266  0.32354  1.8818
+

While the other parameters converge to credible values, the variance +of k2 (omega2.k2) converges to a very small value. The +printout of the saem.mmkin model shows that the estimated +standard deviation of k2 across the population of soils +(SD.k2) is ill-defined, indicating overparameterisation of +this model.

+

When the DFOP model is fitted with the two-component error model, we +also observe that the estimated variance of k2 becomes very small, while +being ill-defined, as illustrated by the excessive confidence interval +of SD.k2.

+
+f_parent_saemix_dfop_tc <- mkin::saem(f_parent_mkin_tc["DFOP", ], quiet = TRUE,
+  control = saemix_control, transformations = "saemix")
+f_parent_saemix_dfop_tc_moreiter <- mkin::saem(f_parent_mkin_tc["DFOP", ], quiet = TRUE,
+  control = saemix_control_moreiter, transformations = "saemix")
+plot(f_parent_saemix_dfop_tc$so, plot.type = "convergence")
+

+
+print(f_parent_saemix_dfop_tc)
+
Kinetic nonlinear mixed-effects model fit by SAEM
+Structural model:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+           time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+           * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Likelihood computed by importance sampling
+  AIC BIC logLik
+  666 664   -323
+
+Fitted parameters:
+          estimate     lower    upper
+DMTA_0    98.24165  96.29190 100.1914
+k1         0.06421   0.03352   0.0949
+k2         0.00866   0.00617   0.0111
+g          0.95340   0.91218   0.9946
+a.1        1.06463   0.87979   1.2495
+b.1        0.02964   0.02266   0.0366
+SD.DMTA_0  2.03611   0.40361   3.6686
+SD.k1      0.59534   0.25692   0.9338
+SD.k2      0.00042 -73.00540  73.0062
+SD.g       1.04234   0.37189   1.7128
+

Doubling the number of iterations in the first phase of the algorithm +leads to a slightly lower likelihood, and therefore to slightly higher +AIC and BIC values. With even more iterations, the algorithm stops with +an error message. This is related to the variance of k2 approximating +zero and has been submitted as a bug +to the saemix package, as the algorithm does not converge in this +case.

+

An alternative way to fit DFOP in combination with the two-component +error model is to use the model formulation with transformed parameters +as used per default in mkin. When using this option, convergence is +slower, but eventually the algorithm stops as well with the same error +message.

+

The four combinations (SFO/const, SFO/tc, DFOP/const and DFOP/tc) and +the version with increased iterations can be compared using the model +comparison function of the saemix package:

+
+AIC_parent_saemix <- saemix::compare.saemix(
+  f_parent_saemix_sfo_const$so,
+  f_parent_saemix_sfo_tc$so,
+  f_parent_saemix_dfop_const$so,
+  f_parent_saemix_dfop_tc$so,
+  f_parent_saemix_dfop_tc_moreiter$so)
+
Likelihoods calculated by importance sampling
+
+rownames(AIC_parent_saemix) <- c(
+  "SFO const", "SFO tc", "DFOP const", "DFOP tc", "DFOP tc more iterations")
+print(AIC_parent_saemix)
+
                           AIC    BIC
+SFO const               796.38 795.34
+SFO tc                  798.38 797.13
+DFOP const              705.75 703.88
+DFOP tc                 665.67 663.59
+DFOP tc more iterations 665.85 663.76
+

In order to check the influence of the likelihood calculation +algorithms implemented in saemix, the likelihood from Gaussian +quadrature is added to the best fit, and the AIC values obtained from +the three methods are compared.

+
+f_parent_saemix_dfop_tc$so <-
+  saemix::llgq.saemix(f_parent_saemix_dfop_tc$so)
+AIC_parent_saemix_methods <- c(
+  is = AIC(f_parent_saemix_dfop_tc$so, method = "is"),
+  gq = AIC(f_parent_saemix_dfop_tc$so, method = "gq"),
+  lin = AIC(f_parent_saemix_dfop_tc$so, method = "lin")
+)
+print(AIC_parent_saemix_methods)
+
    is     gq    lin 
+665.67 665.74 665.13
+

The AIC values based on importance sampling and Gaussian quadrature +are very similar. Using linearisation is known to be less accurate, but +still gives a similar value.

+

In order to illustrate that the comparison of the three method +depends on the degree of convergence obtained in the fit, the same +comparison is shown below for the fit using the defaults for the number +of iterations and the number of MCMC chains.

+

When using OpenBlas for linear algebra, there is a large difference +in the values obtained with Gaussian quadrature, so the larger number of +iterations makes a lot of difference. When using the LAPACK version +coming with Debian Bullseye, the AIC based on Gaussian quadrature is +almost the same as the one obtained with the other methods, also when +using defaults for the fit.

+
+f_parent_saemix_dfop_tc_defaults <- mkin::saem(f_parent_mkin_tc["DFOP", ])
+f_parent_saemix_dfop_tc_defaults$so <-
+  saemix::llgq.saemix(f_parent_saemix_dfop_tc_defaults$so)
+AIC_parent_saemix_methods_defaults <- c(
+  is = AIC(f_parent_saemix_dfop_tc_defaults$so, method = "is"),
+  gq = AIC(f_parent_saemix_dfop_tc_defaults$so, method = "gq"),
+  lin = AIC(f_parent_saemix_dfop_tc_defaults$so, method = "lin")
+)
+print(AIC_parent_saemix_methods_defaults)
+
    is     gq    lin 
+670.09 669.37 671.29
+
+
+
+

Comparison +

+

The following table gives the AIC values obtained with both backend +packages using the same control parameters (800 iterations burn-in, 300 +iterations second phase, 15 chains).

+
+AIC_all <- data.frame(
+  check.names = FALSE,
+  "Degradation model" = c("SFO", "SFO", "DFOP", "DFOP"),
+  "Error model" = c("const", "tc", "const", "tc"),
+  nlme = c(AIC(f_parent_nlme_sfo_const), AIC(f_parent_nlme_sfo_tc), NA, AIC(f_parent_nlme_dfop_tc)),
+  saemix_lin = sapply(list(f_parent_saemix_sfo_const$so, f_parent_saemix_sfo_tc$so,
+    f_parent_saemix_dfop_const$so, f_parent_saemix_dfop_tc$so), AIC, method = "lin"),
+  saemix_is = sapply(list(f_parent_saemix_sfo_const$so, f_parent_saemix_sfo_tc$so,
+    f_parent_saemix_dfop_const$so, f_parent_saemix_dfop_tc$so), AIC, method = "is")
+)
+kable(AIC_all)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Degradation modelError modelnlmesaemix_linsaemix_is
SFOconst796.60796.60796.38
SFOtc798.60798.60798.38
DFOPconstNA704.95705.75
DFOPtc671.91665.13665.67
+
+
+
+

Conclusion +

+

A more detailed analysis of the dimethenamid dataset confirmed that +the DFOP model provides the most appropriate description of the decline +of the parent compound in these data. On the other hand, closer +inspection of the results revealed that the variability of the k2 +parameter across the population of soils is ill-defined. This coincides +with the observation that this parameter cannot robustly be quantified +for some of the soils.

+

Regarding the regulatory use of these data, it is claimed that an +improved characterisation of the mean parameter values across the +population is obtained using the nonlinear mixed-effects models +presented here. However, attempts to quantify the variability of the +slower rate constant of the biphasic decline of dimethenamid indicate +that the data are not sufficient to characterise this variability to a +satisfactory precision.

+
+
+

Session Info +

+ +
R version 4.4.2 (2024-10-31)
+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
+
+locale:
+ [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
+ [3] LC_TIME=C                  LC_COLLATE=de_DE.UTF-8    
+ [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
+ [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
+ [9] LC_ADDRESS=C               LC_TELEPHONE=C            
+[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
+
+time zone: Europe/Berlin
+tzcode source: system (glibc)
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+[1] saemix_3.3   npde_3.5     nlme_3.1-166 mkin_1.2.10  knitr_1.49  
+
+loaded via a namespace (and not attached):
+ [1] gtable_0.3.6      jsonlite_1.8.9    dplyr_1.1.4       compiler_4.4.2   
+ [5] tidyselect_1.2.1  parallel_4.4.2    gridExtra_2.3     jquerylib_0.1.4  
+ [9] systemfonts_1.1.0 scales_1.3.0      textshaping_0.4.1 yaml_2.3.10      
+[13] fastmap_1.2.0     lattice_0.22-6    ggplot2_3.5.1     R6_2.5.1         
+[17] generics_0.1.3    lmtest_0.9-40     MASS_7.3-61       htmlwidgets_1.6.4
+[21] tibble_3.2.1      desc_1.4.3        munsell_0.5.1     bslib_0.8.0      
+[25] pillar_1.9.0      rlang_1.1.4       utf8_1.2.4        cachem_1.1.0     
+[29] xfun_0.49         fs_1.6.5          sass_0.4.9        cli_3.6.3        
+[33] pkgdown_2.1.1     magrittr_2.0.3    digest_0.6.37     grid_4.4.2       
+[37] mclust_6.1.1      lifecycle_1.0.4   vctrs_0.6.5       evaluate_1.0.1   
+[41] glue_1.8.0        codetools_0.2-20  ragg_1.3.3        zoo_1.8-12       
+[45] fansi_1.0.6       colorspace_2.1-1  rmarkdown_2.29    pkgconfig_2.0.3  
+[49] tools_4.4.2       htmltools_0.5.8.1
+
+
+

References +

+ +
+
+EFSA. 2018. “Peer Review of the Pesticide Risk Assessment of the +Active Substance Dimethenamid-p.” EFSA Journal 16: 5211. +
+
+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. +
+
+Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria. +2018. “Renewal Assessment Report +Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour, Rev. 2 - +November 2017.” https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716. +
+
+
+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const-1.png new file mode 100644 index 00000000..627e5c95 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const_test-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const_test-1.png new file mode 100644 index 00000000..627e5c95 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const_test-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_tc_test-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_tc_test-1.png new file mode 100644 index 00000000..9f40fc35 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_tc_test-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_sfo_const-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_sfo_const-1.png new file mode 100644 index 00000000..0ed7448d Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_mkin_sfo_const-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_const-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_const-1.png new file mode 100644 index 00000000..84a6fc92 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_const-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_tc-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_tc-1.png new file mode 100644 index 00000000..fa5d34f0 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_tc-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_const-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_const-1.png new file mode 100644 index 00000000..7862fc65 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_const-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_tc-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_tc-1.png new file mode 100644 index 00000000..d941f3e6 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_tc-1.png differ diff --git a/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/plot_parent_nlme-1.png b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/plot_parent_nlme-1.png new file mode 100644 index 00000000..6bcf3434 Binary files /dev/null and b/docs/dev/articles/web_only/dimethenamid_2018_files/figure-html/plot_parent_nlme-1.png differ diff --git a/docs/dev/articles/web_only/multistart.html b/docs/dev/articles/web_only/multistart.html new file mode 100644 index 00000000..83073979 --- /dev/null +++ b/docs/dev/articles/web_only/multistart.html @@ -0,0 +1,179 @@ + + + + + + + +Short demo of the multistart method • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Short demo of the multistart method

+

Johannes +Ranke

+ +

Last change 20 April 2023 +(rebuilt 2025-02-14)

+ + Source: vignettes/web_only/multistart.rmd +
multistart.rmd
+
+ + + +

The dimethenamid data from 2018 from seven soils is used as example +data in this vignette.

+
+library(mkin)
+dmta_ds <- lapply(1:7, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+  ds_i
+})
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+dmta_ds[["Elliot 1"]] <- dmta_ds[["Elliot 2"]] <- NULL
+

First, we check the DFOP model with the two-component error model and +random effects for all degradation parameters.

+
+f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE)
+f_saem_full <- saem(f_mmkin)
+illparms(f_saem_full)
+
## [1] "sd(log_k2)"
+

We see that not all variability parameters are identifiable. The +illparms function tells us that the confidence interval for +the standard deviation of ‘log_k2’ includes zero. We check this +assessment using multiple runs with different starting values.

+
+f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16)
+parplot(f_saem_full_multi, lpos = "topleft")
+

+

This confirms that the variance of k2 is the most problematic +parameter, so we reduce the parameter distribution model by removing the +intersoil variability for k2.

+
+f_saem_reduced <- stats::update(f_saem_full, no_random_effect = "log_k2")
+illparms(f_saem_reduced)
+f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16)
+parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2))
+

+

The results confirm that all remaining parameters can be determined +with sufficient certainty.

+

We can also analyse the log-likelihoods obtained in the multiple +runs:

+
+llhist(f_saem_reduced_multi)
+

+

We can use the anova method to compare the models.

+
+anova(f_saem_full, best(f_saem_full_multi),
+  f_saem_reduced, best(f_saem_reduced_multi), test = TRUE)
+
## Data: 155 observations of 1 variable(s) grouped in 6 datasets
+## 
+##                            npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)
+## f_saem_reduced                9 663.67 661.80 -322.84                     
+## best(f_saem_reduced_multi)    9 663.65 661.78 -322.82 0.0219  0           
+## f_saem_full                  10 670.09 668.01 -325.05 0.0000  1          1
+## best(f_saem_full_multi)      10 665.61 663.52 -322.80 4.4870  0
+

The reduced model results in lower AIC and BIC values, so it is +clearly preferable. Using multiple starting values gives a large +improvement in case of the full model, because it is less well-defined, +which impedes convergence. For the reduced model, using multiple +starting values only results in a small improvement of the model +fit.

+
+
+ + + +
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 2.1.1.

+
+ +
+
+ + + + + + + diff --git a/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png new file mode 100644 index 00000000..19b68cfe Binary files /dev/null and b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-4-1.png b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-4-1.png new file mode 100644 index 00000000..034b170c Binary files /dev/null and b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-4-1.png differ diff --git a/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-5-1.png b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-5-1.png new file mode 100644 index 00000000..c8e918cd Binary files /dev/null and b/docs/dev/articles/web_only/multistart_files/figure-html/unnamed-chunk-5-1.png differ diff --git a/docs/dev/articles/web_only/saem_benchmarks.html b/docs/dev/articles/web_only/saem_benchmarks.html new file mode 100644 index 00000000..d4cd0d6f --- /dev/null +++ b/docs/dev/articles/web_only/saem_benchmarks.html @@ -0,0 +1,789 @@ + + + + + + + +Benchmark timings for saem.mmkin • mkin + + + + + + + + + Skip to contents + + +
+ + + + +
+
+ +

Benchmark timings for saem.mmkin

+

Johannes +Ranke

+ +

Last change 17 February 2023 +(rebuilt 2025-02-14)

+ + Source: vignettes/web_only/saem_benchmarks.rmd +
saem_benchmarks.rmd
+
+ + + +

Each system is characterized by operating system type, CPU type, mkin +version, saemix version and R version. A compiler was available, so if +no analytical solution was available, compiled ODE models are used.

+

Every fit is only performed once, so the accuracy of the benchmarks +is limited.

+

For the initial mmkin fits, we use all available cores.

+
+n_cores <- parallel::detectCores()
+
+

Test data +

+

Please refer to the vignette dimethenamid_2018 for an +explanation of the following preprocessing.

+
+dmta_ds <- lapply(1:7, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+  ds_i
+})
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL
+
+
+

Test cases +

+
+

Parent only +

+
+parent_mods <- c("SFO", "DFOP", "SFORB", "HS")
+parent_sep_const <- mmkin(parent_mods, dmta_ds, quiet = TRUE, cores = n_cores)
+parent_sep_tc <- update(parent_sep_const, error_model = "tc")
+
+t1 <- system.time(sfo_const <- saem(parent_sep_const["SFO", ]))[["elapsed"]]
+t2 <- system.time(dfop_const <- saem(parent_sep_const["DFOP", ]))[["elapsed"]]
+t3 <- system.time(sforb_const <- saem(parent_sep_const["SFORB", ]))[["elapsed"]]
+t4 <- system.time(hs_const <- saem(parent_sep_const["HS", ]))[["elapsed"]]
+t5 <- system.time(sfo_tc <- saem(parent_sep_tc["SFO", ]))[["elapsed"]]
+t6 <- system.time(dfop_tc <- saem(parent_sep_tc["DFOP", ]))[["elapsed"]]
+t7 <- system.time(sforb_tc <- saem(parent_sep_tc["SFORB", ]))[["elapsed"]]
+t8 <- system.time(hs_tc <- saem(parent_sep_tc["HS", ]))[["elapsed"]]
+
+anova(
+  sfo_const, dfop_const, sforb_const, hs_const,
+  sfo_tc, dfop_tc, sforb_tc, hs_tc) |> kable(, digits = 1)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
nparAICBICLik
sfo_const5796.3795.3-393.2
sfo_tc6798.3797.1-393.2
dfop_const9709.4707.5-345.7
sforb_const9710.0708.1-346.0
hs_const9713.7711.8-347.8
dfop_tc10670.1668.0-325.0
sforb_tc10662.9660.8-321.4
hs_tc10667.2665.1-323.6
+

The above model comparison suggests to use the SFORB model with +two-component error. For comparison, we keep the DFOP model with +two-component error, as it competes with SFORB for biphasic curves.

+
+illparms(dfop_tc)
+
## [1] "sd(log_k2)"
+
+illparms(sforb_tc)
+
## [1] "sd(log_k_DMTA_bound_free)"
+

For these two models, random effects for the transformed parameters +k2 and k_DMTA_bound_free could not be +quantified.

+
+
+

One metabolite +

+

We remove parameters that were found to be ill-defined in the parent +only fits.

+
+one_met_mods <- list(
+  DFOP_SFO = mkinmod(
+    DMTA = mkinsub("DFOP", "M23"),
+    M23 = mkinsub("SFO")),
+  SFORB_SFO = mkinmod(
+    DMTA = mkinsub("SFORB", "M23"),
+    M23 = mkinsub("SFO")))
+
+one_met_sep_const <- mmkin(one_met_mods, dmta_ds, error_model = "const",
+  cores = n_cores, quiet = TRUE)
+one_met_sep_tc <- mmkin(one_met_mods, dmta_ds, error_model = "tc",
+  cores = n_cores, quiet = TRUE)
+
+t9 <- system.time(dfop_sfo_tc <- saem(one_met_sep_tc["DFOP_SFO", ],
+    no_random_effect = "log_k2"))[["elapsed"]]
+t10 <- system.time(sforb_sfo_tc <- saem(one_met_sep_tc["SFORB_SFO", ],
+    no_random_effect = "log_k_DMTA_bound_free"))[["elapsed"]]
+
+
+

Three metabolites +

+

For the case of three metabolites, we only keep the SFORB model in +order to limit the time for compiling this vignette, and as fitting in +parallel may disturb the benchmark. Again, we do not include random +effects that were ill-defined in previous fits of subsets of the +degradation model.

+
+illparms(sforb_sfo_tc)
+
+three_met_mods <- list(
+  SFORB_SFO3_plus = mkinmod(
+    DMTA = mkinsub("SFORB", c("M23", "M27", "M31")),
+    M23 = mkinsub("SFO"),
+    M27 = mkinsub("SFO"),
+    M31 = mkinsub("SFO", "M27", sink = FALSE)))
+
+three_met_sep_tc <- mmkin(three_met_mods, dmta_ds, error_model = "tc",
+  cores = n_cores, quiet = TRUE)
+
+t11 <- system.time(sforb_sfo3_plus_const <- saem(three_met_sep_tc["SFORB_SFO3_plus", ],
+    no_random_effect = "log_k_DMTA_bound_free"))[["elapsed"]]
+
+
+
+

Results +

+

Benchmarks for all available error models are shown. They are +intended for improving mkin, not for comparing CPUs or operating +systems. All trademarks belong to their respective owners.

+
+

Parent only +

+

Constant variance for SFO, DFOP, SFORB and HS.

+ ++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
CPUOSmkinsaemixt1t2t3t4
Ryzen 7 1700Linux1.2.03.22.1404.6264.3284.998
Ryzen 7 1700Linux1.2.23.22.4274.5504.2174.851
Ryzen 9 7950XLinux1.2.13.21.3522.8132.4012.074
Ryzen 9 7950XLinux1.2.23.21.3282.7382.3362.023
Ryzen 9 7950XLinux1.2.33.21.1182.0362.0102.088
Ryzen 9 7950XLinux1.2.33.21.4192.3741.9262.398
Ryzen 9 7950XLinux1.2.43.20.9722.5501.9872.055
Intel(R) Xeon(R) Gold 6134 CPU @ 3.20GHzLinux1.2.63.22.9986.5236.1264.721
Ryzen 9 7950XLinux1.2.63.21.1352.0252.4062.478
Ryzen 9 7950XLinux1.2.93.31.0861.9911.9492.411
Ryzen 9 7950XLinux1.2.103.31.1152.2771.9452.134
+

Two-component error fits for SFO, DFOP, SFORB and HS.

+ ++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
CPUOSmkinsaemixt5t6t7t8
Ryzen 7 1700Linux1.2.03.25.6787.4418.0007.980
Ryzen 7 1700Linux1.2.23.25.3527.2018.1748.401
Ryzen 9 7950XLinux1.2.13.22.3883.0333.5323.310
Ryzen 9 7950XLinux1.2.23.22.3412.9683.4653.341
Ryzen 9 7950XLinux1.2.33.22.1593.5843.3073.460
Ryzen 9 7950XLinux1.2.33.22.3483.1343.2533.530
Ryzen 9 7950XLinux1.2.43.22.1273.5873.4333.595
Intel(R) Xeon(R) Gold 6134 CPU @ 3.20GHzLinux1.2.63.25.0708.4648.5257.599
Ryzen 9 7950XLinux1.2.63.22.1613.3253.6693.153
Ryzen 9 7950XLinux1.2.93.32.4263.1963.2563.322
Ryzen 9 7950XLinux1.2.103.32.3723.1373.1003.281
+
+
+

One metabolite +

+

Two-component error for DFOP-SFO and SFORB-SFO.

+ ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
CPUOSmkinsaemixt9t10
Ryzen 7 1700Linux1.2.03.224.465800.266
Ryzen 7 1700Linux1.2.23.225.193798.580
Ryzen 9 7950XLinux1.2.13.211.247285.216
Ryzen 9 7950XLinux1.2.23.211.242284.258
Ryzen 9 7950XLinux1.2.33.211.796216.012
Ryzen 9 7950XLinux1.2.33.212.841292.688
Ryzen 9 7950XLinux1.2.43.212.160265.934
Intel(R) Xeon(R) Gold 6134 CPU @ 3.20GHzLinux1.2.63.230.168748.675
Ryzen 9 7950XLinux1.2.63.212.007286.757
Ryzen 9 7950XLinux1.2.93.312.420289.338
Ryzen 9 7950XLinux1.2.103.311.590279.782
+
+
+

Three metabolites +

+

Two-component error for SFORB-SFO3-plus

+ +++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
CPUOSmkinsaemixt11
Ryzen 7 1700Linux1.2.03.21289.198
Ryzen 7 1700Linux1.2.23.21312.445
Ryzen 9 7950XLinux1.2.13.2489.939
Ryzen 9 7950XLinux1.2.23.2482.970
Ryzen 9 7950XLinux1.2.33.2392.364
Ryzen 9 7950XLinux1.2.33.2483.027
Ryzen 9 7950XLinux1.2.43.2456.252
Intel(R) Xeon(R) Gold 6134 CPU @ 3.20GHzLinux1.2.63.21235.028
Ryzen 9 7950XLinux1.2.63.2480.577
Ryzen 9 7950XLinux1.2.93.3485.836
Ryzen 9 7950XLinux1.2.103.3469.036
+
+
+
+
+ + + + +
+ + + + + + + -- cgit v1.2.1