From 218a9c55bd80fb708b15fa7196422f759bfe4b27 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 May 2020 16:20:23 +0200 Subject: Further formatting improvement of benchmark vignette Also, use .rmd extension instead of .Rmd for vignettes. --- vignettes/web_only/FOCUS_Z.Rmd | 255 --------------------------------- vignettes/web_only/FOCUS_Z.rmd | 255 +++++++++++++++++++++++++++++++++ vignettes/web_only/NAFTA_examples.Rmd | 243 ------------------------------- vignettes/web_only/NAFTA_examples.rmd | 243 +++++++++++++++++++++++++++++++ vignettes/web_only/benchmarks.Rmd | 150 ------------------- vignettes/web_only/benchmarks.html | 136 ++++++++---------- vignettes/web_only/benchmarks.rmd | 157 ++++++++++++++++++++ vignettes/web_only/compiled_models.Rmd | 141 ------------------ vignettes/web_only/compiled_models.rmd | 141 ++++++++++++++++++ vignettes/web_only/mkin_benchmarks.rda | Bin 854 -> 854 bytes 10 files changed, 852 insertions(+), 869 deletions(-) delete mode 100644 vignettes/web_only/FOCUS_Z.Rmd create mode 100644 vignettes/web_only/FOCUS_Z.rmd delete mode 100644 vignettes/web_only/NAFTA_examples.Rmd create mode 100644 vignettes/web_only/NAFTA_examples.rmd delete mode 100644 vignettes/web_only/benchmarks.Rmd create mode 100644 vignettes/web_only/benchmarks.rmd delete mode 100644 vignettes/web_only/compiled_models.Rmd create mode 100644 vignettes/web_only/compiled_models.rmd (limited to 'vignettes/web_only') diff --git a/vignettes/web_only/FOCUS_Z.Rmd b/vignettes/web_only/FOCUS_Z.Rmd deleted file mode 100644 index 2da7fde7..00000000 --- a/vignettes/web_only/FOCUS_Z.Rmd +++ /dev/null @@ -1,255 +0,0 @@ ---- -title: Example evaluation of FOCUS dataset Z -author: Johannes Ranke -date: "`r Sys.Date()`" -output: - html_document: - toc: true - toc_float: true - code_folding: show - fig_retina: null -bibliography: ../references.bib -vignette: > - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -[Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany](http://www.jrwb.de)
-[Privatdozent at the University of Bremen](http://chem.uft.uni-bremen.de/ranke) - -```{r, include = FALSE} -require(knitr) -options(digits = 5) -opts_chunk$set(engine='R', tidy = FALSE) -``` - -# The data - -The following code defines the example dataset from Appendix 7 to the FOCUS kinetics -report [@FOCUSkinetics2014, p. 354]. - -```{r, echo = TRUE, fig = TRUE, fig.width = 8, fig.height = 7} -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). - -```{r FOCUS_2006_Z_fits_1, echo=TRUE, fig.height=6} -Z.2a <- mkinmod(Z0 = mkinsub("SFO", "Z1"), - Z1 = mkinsub("SFO")) -m.Z.2a <- mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE) -plot_sep(m.Z.2a) -summary(m.Z.2a, data = FALSE)$bpar -``` - -As obvious from the parameter summary (the \texttt{bpar} 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: - -```{r FOCUS_2006_Z_fits_2, echo=TRUE, fig.height=6} -Z.2a.ff <- mkinmod(Z0 = mkinsub("SFO", "Z1"), - Z1 = mkinsub("SFO"), - use_of_ff = "max") - -m.Z.2a.ff <- mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE) -plot_sep(m.Z.2a.ff) -summary(m.Z.2a.ff, data = FALSE)$bpar -``` - -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. -\footnote{If the model formulation without formation fractions -is used, the same effect can be obtained by fixing the parameter \texttt{k\_Z\_sink} -to a value of zero.} - -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. - -```{r FOCUS_2006_Z_fits_3, echo=TRUE, fig.height=6} -Z.3 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), - Z1 = mkinsub("SFO"), use_of_ff = "max") -m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE) -plot_sep(m.Z.3) -summary(m.Z.3, data = FALSE)$bpar -``` - -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. - -```{r FOCUS_2006_Z_fits_5, echo=TRUE, fig.height=7} -Z.5 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), - Z1 = mkinsub("SFO", "Z2", sink = FALSE), - Z2 = mkinsub("SFO"), use_of_ff = "max") -m.Z.5 <- mkinfit(Z.5, FOCUS_2006_Z_mkin, quiet = TRUE) -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. - -```{r FOCUS_2006_Z_fits_6, echo=TRUE, fig.height=8} -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") -m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, - parms.ini = m.Z.5$bparms.ode, - quiet = TRUE) -plot_sep(m.Z.FOCUS) -summary(m.Z.FOCUS, data = FALSE)$bpar -endpoints(m.Z.FOCUS) -``` - -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 $\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. - -```{r FOCUS_2006_Z_fits_7, echo=TRUE, fig.height=8} -Z.mkin.1 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), - Z1 = mkinsub("SFO", "Z2", sink = FALSE), - Z2 = mkinsub("SFO", "Z3"), - Z3 = mkinsub("SFORB")) -m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE) -plot_sep(m.Z.mkin.1) -summary(m.Z.mkin.1, data = FALSE)$cov.unscaled -``` - -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. - -```{r FOCUS_2006_Z_fits_9, echo=TRUE, fig.height=8} -Z.mkin.3 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), - Z1 = mkinsub("SFO", "Z2", sink = FALSE), - Z2 = mkinsub("SFO")) -m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE) -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. - -```{r FOCUS_2006_Z_fits_10, echo=TRUE, fig.height=8} -Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), - Z1 = mkinsub("SFO", "Z2", sink = FALSE), - Z2 = mkinsub("SFO", "Z3"), - Z3 = mkinsub("SFO")) -m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, - parms.ini = m.Z.mkin.3$bparms.ode, - quiet = TRUE) -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. - -```{r FOCUS_2006_Z_fits_11, echo=TRUE, fig.height=8} -Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), - Z1 = mkinsub("SFO", "Z2", sink = FALSE), - Z2 = mkinsub("SFO", "Z3"), - Z3 = mkinsub("SFORB")) -m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, - parms.ini = m.Z.mkin.4$bparms.ode[1:4], - quiet = TRUE) -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 -\texttt{k\_Z3\_bound\_free} is excessively small, it seems reasonable to fix it to -zero. - -```{r FOCUS_2006_Z_fits_11a, echo=TRUE} -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) -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 -\texttt{Z.mkin.5a} 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. - -```{r FOCUS_2006_Z_fits_11b, echo=TRUE} -mkinparplot(m.Z.mkin.5a) -``` - -The endpoints obtained with this model are - -```{r FOCUS_2006_Z_fits_11b_endpoints, echo=TRUE} -endpoints(m.Z.mkin.5a) -``` - -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 - - diff --git a/vignettes/web_only/FOCUS_Z.rmd b/vignettes/web_only/FOCUS_Z.rmd new file mode 100644 index 00000000..2da7fde7 --- /dev/null +++ b/vignettes/web_only/FOCUS_Z.rmd @@ -0,0 +1,255 @@ +--- +title: Example evaluation of FOCUS dataset Z +author: Johannes Ranke +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: true + code_folding: show + fig_retina: null +bibliography: ../references.bib +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +[Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany](http://www.jrwb.de)
+[Privatdozent at the University of Bremen](http://chem.uft.uni-bremen.de/ranke) + +```{r, include = FALSE} +require(knitr) +options(digits = 5) +opts_chunk$set(engine='R', tidy = FALSE) +``` + +# The data + +The following code defines the example dataset from Appendix 7 to the FOCUS kinetics +report [@FOCUSkinetics2014, p. 354]. + +```{r, echo = TRUE, fig = TRUE, fig.width = 8, fig.height = 7} +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). + +```{r FOCUS_2006_Z_fits_1, echo=TRUE, fig.height=6} +Z.2a <- mkinmod(Z0 = mkinsub("SFO", "Z1"), + Z1 = mkinsub("SFO")) +m.Z.2a <- mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE) +plot_sep(m.Z.2a) +summary(m.Z.2a, data = FALSE)$bpar +``` + +As obvious from the parameter summary (the \texttt{bpar} 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: + +```{r FOCUS_2006_Z_fits_2, echo=TRUE, fig.height=6} +Z.2a.ff <- mkinmod(Z0 = mkinsub("SFO", "Z1"), + Z1 = mkinsub("SFO"), + use_of_ff = "max") + +m.Z.2a.ff <- mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE) +plot_sep(m.Z.2a.ff) +summary(m.Z.2a.ff, data = FALSE)$bpar +``` + +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. +\footnote{If the model formulation without formation fractions +is used, the same effect can be obtained by fixing the parameter \texttt{k\_Z\_sink} +to a value of zero.} + +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. + +```{r FOCUS_2006_Z_fits_3, echo=TRUE, fig.height=6} +Z.3 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), + Z1 = mkinsub("SFO"), use_of_ff = "max") +m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE) +plot_sep(m.Z.3) +summary(m.Z.3, data = FALSE)$bpar +``` + +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. + +```{r FOCUS_2006_Z_fits_5, echo=TRUE, fig.height=7} +Z.5 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), + Z1 = mkinsub("SFO", "Z2", sink = FALSE), + Z2 = mkinsub("SFO"), use_of_ff = "max") +m.Z.5 <- mkinfit(Z.5, FOCUS_2006_Z_mkin, quiet = TRUE) +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. + +```{r FOCUS_2006_Z_fits_6, echo=TRUE, fig.height=8} +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") +m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, + parms.ini = m.Z.5$bparms.ode, + quiet = TRUE) +plot_sep(m.Z.FOCUS) +summary(m.Z.FOCUS, data = FALSE)$bpar +endpoints(m.Z.FOCUS) +``` + +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 $\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. + +```{r FOCUS_2006_Z_fits_7, echo=TRUE, fig.height=8} +Z.mkin.1 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), + Z1 = mkinsub("SFO", "Z2", sink = FALSE), + Z2 = mkinsub("SFO", "Z3"), + Z3 = mkinsub("SFORB")) +m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE) +plot_sep(m.Z.mkin.1) +summary(m.Z.mkin.1, data = FALSE)$cov.unscaled +``` + +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. + +```{r FOCUS_2006_Z_fits_9, echo=TRUE, fig.height=8} +Z.mkin.3 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), + Z1 = mkinsub("SFO", "Z2", sink = FALSE), + Z2 = mkinsub("SFO")) +m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE) +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. + +```{r FOCUS_2006_Z_fits_10, echo=TRUE, fig.height=8} +Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), + Z1 = mkinsub("SFO", "Z2", sink = FALSE), + Z2 = mkinsub("SFO", "Z3"), + Z3 = mkinsub("SFO")) +m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, + parms.ini = m.Z.mkin.3$bparms.ode, + quiet = TRUE) +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. + +```{r FOCUS_2006_Z_fits_11, echo=TRUE, fig.height=8} +Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), + Z1 = mkinsub("SFO", "Z2", sink = FALSE), + Z2 = mkinsub("SFO", "Z3"), + Z3 = mkinsub("SFORB")) +m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, + parms.ini = m.Z.mkin.4$bparms.ode[1:4], + quiet = TRUE) +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 +\texttt{k\_Z3\_bound\_free} is excessively small, it seems reasonable to fix it to +zero. + +```{r FOCUS_2006_Z_fits_11a, echo=TRUE} +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) +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 +\texttt{Z.mkin.5a} 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. + +```{r FOCUS_2006_Z_fits_11b, echo=TRUE} +mkinparplot(m.Z.mkin.5a) +``` + +The endpoints obtained with this model are + +```{r FOCUS_2006_Z_fits_11b_endpoints, echo=TRUE} +endpoints(m.Z.mkin.5a) +``` + +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 + + diff --git a/vignettes/web_only/NAFTA_examples.Rmd b/vignettes/web_only/NAFTA_examples.Rmd deleted file mode 100644 index 26a9240a..00000000 --- a/vignettes/web_only/NAFTA_examples.Rmd +++ /dev/null @@ -1,243 +0,0 @@ ---- -title: "Evaluation of example datasets from Attachment 1 to the US EPA SOP for the NAFTA guidance" -author: "Johannes Ranke" -date: "`r Sys.Date()`" -output: - html_document: - toc: true - toc_float: - collapsed: false - mathjax: null - fig_retina: null -references: -- id: usepa2015 - title: Standard Operating Procedure for Using the NAFTA - Guidance to Calculate Representative Half-life Values and Characterizing - Pesticide Degradation - author: - - family: US EPA - type: report - issued: - year: 2015 - url: https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance -vignette: > - %\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -# Introduction - -In this document, the example evaluations provided in Attachment 1 to the SOP -of US EPA for using the NAFTA guidance [@usepa2015] 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. - -```{r, include = FALSE} -library(knitr) -opts_chunk$set(tidy = FALSE, cache = FALSE, fig.height = 7) -library("mkin", quietly = TRUE) -``` - -# 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 - -```{r p5a} -p5a <- nafta(NAFTA_SOP_Attachment[["p5a"]]) -plot(p5a) -print(p5a) -``` - -## Example on page 5, lower panel - -```{r p5b} -p5b <- nafta(NAFTA_SOP_Attachment[["p5b"]]) -plot(p5b) -print(p5b) -``` - -## Example on page 6 - -```{r p6} -p6 <- nafta(NAFTA_SOP_Attachment[["p6"]]) -plot(p6) -print(p6) -``` - -## Example on page 7 - -```{r p7} -p7 <- nafta(NAFTA_SOP_Attachment[["p7"]]) -plot(p7) -print(p7) -``` - -# 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. - -```{r p8} -p8 <- nafta(NAFTA_SOP_Attachment[["p8"]], parms.ini = c(k__iore_parent_sink = 1e-3)) -plot(p8) -print(p8) -``` - -# Examples where SFO was not selected for an abiotic study - -## Example on page 9, upper panel - -```{r p9a} -p9a <- nafta(NAFTA_SOP_Attachment[["p9a"]]) -plot(p9a) -print(p9a) -``` - -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 - -```{r p9b} -p9b <- nafta(NAFTA_SOP_Attachment[["p9b"]]) -plot(p9b) -print(p9b) -``` - -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 - -```{r p10} -p10 <- nafta(NAFTA_SOP_Attachment[["p10"]]) -plot(p10) -print(p10) -``` - -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 - -```{r p11} -p11 <- nafta(NAFTA_SOP_Attachment[["p11"]]) -plot(p11) -print(p11) -``` - -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 - -```{r p12a} -p12a <- nafta(NAFTA_SOP_Attachment[["p12a"]]) -plot(p12a) -print(p12a) -``` - -## Example on page 12, lower panel - -```{r p12b} -p12b <- nafta(NAFTA_SOP_Attachment[["p12b"]]) -plot(p12b) -print(p12b) -``` - -## Example on page 13 - -```{r p13} -p13 <- nafta(NAFTA_SOP_Attachment[["p13"]]) -plot(p13) -print(p13) -``` - -# DT50 not observed in the study and DFOP problems in PestDF - -```{r p14} -p14 <- nafta(NAFTA_SOP_Attachment[["p14"]]) -plot(p14) -print(p14) -``` - -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 - -```{r p15a} -p15a <- nafta(NAFTA_SOP_Attachment[["p15a"]]) -plot(p15a) -print(p15a) -``` - -```{r p15b} -p15b <- nafta(NAFTA_SOP_Attachment[["p15b"]]) -plot(p15b) -print(p15b) -``` - -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 - -```{r p16} -p16 <- nafta(NAFTA_SOP_Attachment[["p16"]]) -plot(p16) -print(p16) -``` - -In PestDF, the DFOP fit seems to have stuck in a local minimum, as mkin finds -a solution with a much lower $\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 diff --git a/vignettes/web_only/NAFTA_examples.rmd b/vignettes/web_only/NAFTA_examples.rmd new file mode 100644 index 00000000..26a9240a --- /dev/null +++ b/vignettes/web_only/NAFTA_examples.rmd @@ -0,0 +1,243 @@ +--- +title: "Evaluation of example datasets from Attachment 1 to the US EPA SOP for the NAFTA guidance" +author: "Johannes Ranke" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: + collapsed: false + mathjax: null + fig_retina: null +references: +- id: usepa2015 + title: Standard Operating Procedure for Using the NAFTA + Guidance to Calculate Representative Half-life Values and Characterizing + Pesticide Degradation + author: + - family: US EPA + type: report + issued: + year: 2015 + url: https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance +vignette: > + %\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +# Introduction + +In this document, the example evaluations provided in Attachment 1 to the SOP +of US EPA for using the NAFTA guidance [@usepa2015] 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. + +```{r, include = FALSE} +library(knitr) +opts_chunk$set(tidy = FALSE, cache = FALSE, fig.height = 7) +library("mkin", quietly = TRUE) +``` + +# 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 + +```{r p5a} +p5a <- nafta(NAFTA_SOP_Attachment[["p5a"]]) +plot(p5a) +print(p5a) +``` + +## Example on page 5, lower panel + +```{r p5b} +p5b <- nafta(NAFTA_SOP_Attachment[["p5b"]]) +plot(p5b) +print(p5b) +``` + +## Example on page 6 + +```{r p6} +p6 <- nafta(NAFTA_SOP_Attachment[["p6"]]) +plot(p6) +print(p6) +``` + +## Example on page 7 + +```{r p7} +p7 <- nafta(NAFTA_SOP_Attachment[["p7"]]) +plot(p7) +print(p7) +``` + +# 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. + +```{r p8} +p8 <- nafta(NAFTA_SOP_Attachment[["p8"]], parms.ini = c(k__iore_parent_sink = 1e-3)) +plot(p8) +print(p8) +``` + +# Examples where SFO was not selected for an abiotic study + +## Example on page 9, upper panel + +```{r p9a} +p9a <- nafta(NAFTA_SOP_Attachment[["p9a"]]) +plot(p9a) +print(p9a) +``` + +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 + +```{r p9b} +p9b <- nafta(NAFTA_SOP_Attachment[["p9b"]]) +plot(p9b) +print(p9b) +``` + +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 + +```{r p10} +p10 <- nafta(NAFTA_SOP_Attachment[["p10"]]) +plot(p10) +print(p10) +``` + +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 + +```{r p11} +p11 <- nafta(NAFTA_SOP_Attachment[["p11"]]) +plot(p11) +print(p11) +``` + +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 + +```{r p12a} +p12a <- nafta(NAFTA_SOP_Attachment[["p12a"]]) +plot(p12a) +print(p12a) +``` + +## Example on page 12, lower panel + +```{r p12b} +p12b <- nafta(NAFTA_SOP_Attachment[["p12b"]]) +plot(p12b) +print(p12b) +``` + +## Example on page 13 + +```{r p13} +p13 <- nafta(NAFTA_SOP_Attachment[["p13"]]) +plot(p13) +print(p13) +``` + +# DT50 not observed in the study and DFOP problems in PestDF + +```{r p14} +p14 <- nafta(NAFTA_SOP_Attachment[["p14"]]) +plot(p14) +print(p14) +``` + +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 + +```{r p15a} +p15a <- nafta(NAFTA_SOP_Attachment[["p15a"]]) +plot(p15a) +print(p15a) +``` + +```{r p15b} +p15b <- nafta(NAFTA_SOP_Attachment[["p15b"]]) +plot(p15b) +print(p15b) +``` + +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 + +```{r p16} +p16 <- nafta(NAFTA_SOP_Attachment[["p16"]]) +plot(p16) +print(p16) +``` + +In PestDF, the DFOP fit seems to have stuck in a local minimum, as mkin finds +a solution with a much lower $\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 diff --git a/vignettes/web_only/benchmarks.Rmd b/vignettes/web_only/benchmarks.Rmd deleted file mode 100644 index 990c2fee..00000000 --- a/vignettes/web_only/benchmarks.Rmd +++ /dev/null @@ -1,150 +0,0 @@ ---- -title: "Benchmark timings for mkin" -author: "Johannes Ranke" -output: - html_document: - toc: true - toc_float: true - code_folding: show - fig_retina: null -date: "`r Sys.Date()`" -vignette: > - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r, include = FALSE} -library(knitr) -opts_chunk$set(tidy = FALSE, cache = FALSE) -library("mkin") -``` - -Each system is characterized by its CPU type, the operating system type and the -mkin version. Currently only values for one system are available. - -```{r} -cpu_model <- benchmarkme::get_cpu()$model_name -operating_system <- Sys.info()[["sysname"]] -mkin_version <- as.character(packageVersion("mkin")) -system_string <- paste0(operating_system, ", ", cpu_model, ", mkin version ", mkin_version) -load("~/git/mkin/vignettes/web_only/mkin_benchmarks.rda") -mkin_benchmarks[system_string, c("CPU", "OS", "mkin")] <- - c(cpu_model, operating_system, mkin_version) - -if (mkin_version > "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: - -```{r parent_only, warning = FALSE} -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: - -```{r 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"), - 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: - -```{r two_metabolites} -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"]] -``` - -```{r results} -mkin_benchmarks[system_string, paste0("t", 1:11)] <- - c(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11) -save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks.rda") -``` - -## Results - -Currently, we only have benchmark information on one system, therefore only the mkin -version is shown with the results below. Timings are in seconds, shorter is better. - -```{r} -rownames(mkin_benchmarks) <- as.character(mkin_benchmarks$mkin) -``` - - -Benchmarks for all available error models are shown. - -### Parent only - -Constant variance and two-component error model: - -```{r} -kable(mkin_benchmarks[, c("t1", "t2")]) -``` - -### One metabolite - -Constant variance, variance by variable and two-component error model: - -```{r} -kable(mkin_benchmarks[, c("t3", "t4", "t5")]) -``` - -### Two metabolites - -Two different datasets, for each constant variance, variance by variable and -two-component error model are shown: - -```{r} -kable(mkin_benchmarks[, paste0("t", 6:11)]) -``` diff --git a/vignettes/web_only/benchmarks.html b/vignettes/web_only/benchmarks.html index 821399e4..043b777f 100644 --- a/vignettes/web_only/benchmarks.html +++ b/vignettes/web_only/benchmarks.html @@ -11,7 +11,7 @@ - + Benchmark timings for mkin @@ -1583,29 +1583,12 @@ div.tocify {

Benchmark timings for mkin

Johannes Ranke

-

2020-05-12

+

2020-05-13

-

Each system is characterized by its CPU type, the operating system type and the mkin version. Currently only values for one system are available.

-
cpu_model <- benchmarkme::get_cpu()$model_name
-operating_system <- Sys.info()[["sysname"]]
-mkin_version <- as.character(packageVersion("mkin"))
-system_string <- paste0(operating_system, ", ", cpu_model, ", mkin version ", mkin_version)
-load("~/git/mkin/vignettes/web_only/mkin_benchmarks.rda")
-mkin_benchmarks[system_string, c("CPU", "OS", "mkin")] <-
-  c(cpu_model, operating_system, mkin_version)
-
-if (mkin_version > "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)
-  }
-}
+

Each system is characterized by its CPU type, the operating system type and the mkin version. Currently only values for one system are available. A compiler was available, so if no analytical solution was available, compiled ODE models are used.

Test cases

Parent only:

@@ -1619,17 +1602,14 @@ t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP

One metabolite:

SFO_SFO <- mkinmod(
   parent = mkinsub("SFO", "m1"),
-  m1 = mkinsub("SFO"))
-
## Successfully compiled differential equation model from auto-generated C code.
-
FOMC_SFO <- mkinmod(
+  m1 = mkinsub("SFO"))
+FOMC_SFO <- mkinmod(
   parent = mkinsub("FOMC", "m1"),
-  m1 = mkinsub("SFO"))
-
## Successfully compiled differential equation model from auto-generated C code.
-
DFOP_SFO <- mkinmod(
+  m1 = mkinsub("SFO"))
+DFOP_SFO <- mkinmod(
   parent = mkinsub("FOMC", "m1"),
-  m1 = mkinsub("SFO"))
-
## Successfully compiled differential equation model from auto-generated C code.
-
t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D)))[["elapsed"]]
+  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),
@@ -1667,127 +1647,123 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks

Results

-

Currently, we only have benchmark information on one system, therefore only the mkin version is shown with the results below. Timings are in seconds, shorter is better.

-
rownames(mkin_benchmarks) <- as.character(mkin_benchmarks$mkin)
+

Currently, we only have benchmark information on one system, therefore only the mkin version is shown with the results below. Timings are in seconds, shorter is better. All results were obtained by serial, i.e. not using multiple computing cores.

Benchmarks for all available error models are shown.

Parent only

-

Constant variance and two-component error model:

-
kable(mkin_benchmarks[, c("t1", "t2")])
+

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

- - - + + + - + - + - + - + - + - - - + + +
t1t2mkin versiont1 [s]t2 [s]
0.9.48.10.9.48.1 3.610 11.019
0.9.49.10.9.49.1 8.184 22.889
0.9.49.20.9.49.2 7.064 12.558
0.9.49.30.9.49.3 7.296 21.239
0.9.49.40.9.49.4 5.936 20.545
0.9.50.21.5593.9290.9.50.21.5473.928

One metabolite

-

Constant variance, variance by variable and two-component error model:

-
kable(mkin_benchmarks[, c("t3", "t4", "t5")])
+

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.

- - - - + + + + - + - + - + - + - + - - - - + + + +
t3t4t5mkin versiont3 [s]t4 [s]t5 [s]
0.9.48.10.9.48.1 3.764 14.347 9.495
0.9.49.10.9.49.1 4.649 13.789 6.395
0.9.49.20.9.49.2 4.786 8.461 5.675
0.9.49.30.9.49.3 4.510 13.805 7.386
0.9.49.40.9.49.4 4.446 15.335 6.002
0.9.50.21.3526.1102.8410.9.50.21.3716.1542.720

Two metabolites

-

Two different datasets, for each constant variance, variance by variable and two-component error model are shown:

-
kable(mkin_benchmarks[, paste0("t", 6:11)])
+

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.

- - - - - - - + + + + + + + - + @@ -1796,7 +1772,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks - + @@ -1805,7 +1781,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks - + @@ -1814,7 +1790,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks - + @@ -1823,7 +1799,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks - + @@ -1832,13 +1808,13 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks - - - - - - - + + + + + + +
t6t7t8t9t10t11mkin versiont6 [s]t7 [s]t8 [s]t9 [s]t10 [s]t11 [s]
0.9.48.10.9.48.1 2.623 4.587 7.52531.267
0.9.49.10.9.49.1 2.542 4.128 4.6325.636
0.9.49.20.9.49.2 2.723 4.478 4.8625.574
0.9.49.30.9.49.3 2.643 4.374 7.0207.365
0.9.49.40.9.49.4 2.635 4.259 4.7375.626
0.9.50.20.7591.2041.2752.8372.0262.9760.9.50.20.7421.2031.2822.8452.0372.987
diff --git a/vignettes/web_only/benchmarks.rmd b/vignettes/web_only/benchmarks.rmd new file mode 100644 index 00000000..7ae12451 --- /dev/null +++ b/vignettes/web_only/benchmarks.rmd @@ -0,0 +1,157 @@ +--- +title: "Benchmark timings for mkin" +author: "Johannes Ranke" +output: + html_document: + toc: true + toc_float: true + code_folding: show + fig_retina: null +date: "`r Sys.Date()`" +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +library(knitr) +opts_chunk$set(tidy = FALSE, cache = FALSE) +library("mkin") +``` + +Each system is characterized by its CPU type, the operating system type and the +mkin version. Currently only values for one system are available. A compiler +was available, so if no analytical solution was available, compiled ODE models are +used. + +```{r include = FALSE} +cpu_model <- benchmarkme::get_cpu()$model_name +operating_system <- Sys.info()[["sysname"]] +mkin_version <- as.character(packageVersion("mkin")) +system_string <- paste0(operating_system, ", ", cpu_model, ", mkin version ", mkin_version) +load("~/git/mkin/vignettes/web_only/mkin_benchmarks.rda") +mkin_benchmarks[system_string, c("CPU", "OS", "mkin")] <- + c(cpu_model, operating_system, mkin_version) + +if (mkin_version > "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: + +```{r parent_only, warning = FALSE} +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: + +```{r one_metabolite, message = FALSE} +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"), + 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: + +```{r two_metabolites, message = FALSE} +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"]] +``` + +```{r results} +mkin_benchmarks[system_string, paste0("t", 1:11)] <- + c(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11) +save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks.rda") +``` + +## Results + +Currently, we only have benchmark information on one system, therefore only the mkin +version is shown with the results below. Timings are in seconds, shorter is better. +All results were obtained by serial, i.e. not using multiple computing cores. + +```{r, include = FALSE} +dimnames(mkin_benchmarks) <- list(as.character(mkin_benchmarks$mkin), + c("CPU", "OS", "mkin version", paste0("t", 1:11, " [s]"))) +``` + + +Benchmarks for all available error models are shown. + +### 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. + +```{r, echo = FALSE} +kable(mkin_benchmarks[, 4:5], rownames = "mkin version") +``` + +### 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. + +```{r, echo = FALSE} +kable(mkin_benchmarks[, 6:8], rownames = "mkin version") +``` + +### 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. + +```{r, echo = FALSE} +kable(mkin_benchmarks[, 9:14], rownames = "mkin version") +``` diff --git a/vignettes/web_only/compiled_models.Rmd b/vignettes/web_only/compiled_models.Rmd deleted file mode 100644 index f99ea808..00000000 --- a/vignettes/web_only/compiled_models.Rmd +++ /dev/null @@ -1,141 +0,0 @@ ---- -title: "Performance benefit by using compiled model definitions in mkin" -author: "Johannes Ranke" -output: - html_document: - toc: true - toc_float: true - code_folding: show - fig_retina: null -date: "`r Sys.Date()`" -vignette: > - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r, include = FALSE} -library(knitr) -opts_chunk$set(tidy = FALSE, cache = FALSE) -``` - -## 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 - -```{r check_gcc, eval = FALSE} -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 - -```{r Rprofile, eval = FALSE} -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 - -```{r HOME, eval = FALSE} -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. - -```{r create_SFO_SFO} -library("mkin", quietly = TRUE) -SFO_SFO <- mkinmod( - parent = mkinsub("SFO", "m1"), - m1 = mkinsub("SFO")) -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. - -```{r benchmark_SFO_SFO, fig.height = 3, message = FALSE, warning = FALSE} -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") -} -``` - -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. - -```{r benchmark_FOMC_SFO, fig.height = 3, warning = FALSE} -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") -} -``` - -Here we get a performance benefit of a factor of -`r factor_FOMC_SFO` -using the version of the differential equation model compiled from C code! - -This vignette was built with mkin `r utils::packageVersion("mkin")` on - -```{r sessionInfo, echo = FALSE} -cat(utils::capture.output(utils::sessionInfo())[1:3], sep = "\n") -if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) { - cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]])) -} -``` diff --git a/vignettes/web_only/compiled_models.rmd b/vignettes/web_only/compiled_models.rmd new file mode 100644 index 00000000..0b8c617a --- /dev/null +++ b/vignettes/web_only/compiled_models.rmd @@ -0,0 +1,141 @@ +--- +title: "Performance benefit by using compiled model definitions in mkin" +author: "Johannes Ranke" +output: + html_document: + toc: true + toc_float: true + code_folding: show + fig_retina: null +date: "`r Sys.Date()`" +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +library(knitr) +opts_chunk$set(tidy = FALSE, cache = FALSE) +``` + +## 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 + +```{r check_gcc, eval = FALSE} +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 + +```{r Rprofile, eval = FALSE} +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 + +```{r HOME, eval = FALSE} +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. + +```{r create_SFO_SFO} +library("mkin", quietly = TRUE) +SFO_SFO <- mkinmod( + parent = mkinsub("SFO", "m1"), + m1 = mkinsub("SFO")) +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. + +```{r benchmark_SFO_SFO, fig.height = 3, message = FALSE, warning = FALSE} +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") +} +``` + +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. + +```{r benchmark_FOMC_SFO, fig.height = 3, warning = FALSE} +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") +} +``` + +Here we get a performance benefit of a factor of +`r factor_FOMC_SFO` +using the version of the differential equation model compiled from C code! + +This vignette was built with mkin `r utils::packageVersion("mkin")` on + +```{r sessionInfo, echo = FALSE} +cat(utils::capture.output(utils::sessionInfo())[1:3], sep = "\n") +if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) { + cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]])) +} +``` diff --git a/vignettes/web_only/mkin_benchmarks.rda b/vignettes/web_only/mkin_benchmarks.rda index 268e5efe..c88ca4d6 100644 Binary files a/vignettes/web_only/mkin_benchmarks.rda and b/vignettes/web_only/mkin_benchmarks.rda differ -- cgit v1.2.1