diff options
Diffstat (limited to 'vignettes')
| -rw-r--r-- | vignettes/FOCUS_D.rmd (renamed from vignettes/FOCUS_D.Rmd) | 146 | ||||
| -rw-r--r-- | vignettes/FOCUS_L.rmd (renamed from vignettes/FOCUS_L.Rmd) | 542 | ||||
| -rw-r--r-- | vignettes/mkin.rmd (renamed from vignettes/mkin.Rmd) | 0 | ||||
| -rw-r--r-- | vignettes/twa.rmd (renamed from vignettes/twa.Rmd) | 0 | ||||
| -rw-r--r-- | vignettes/web_only/FOCUS_Z.rmd (renamed from vignettes/web_only/FOCUS_Z.Rmd) | 0 | ||||
| -rw-r--r-- | vignettes/web_only/NAFTA_examples.rmd (renamed from vignettes/web_only/NAFTA_examples.Rmd) | 0 | ||||
| -rw-r--r-- | vignettes/web_only/benchmarks.html | 136 | ||||
| -rw-r--r-- | vignettes/web_only/benchmarks.rmd (renamed from vignettes/web_only/benchmarks.Rmd) | 39 | ||||
| -rw-r--r-- | vignettes/web_only/compiled_models.rmd (renamed from vignettes/web_only/compiled_models.Rmd) | 282 | ||||
| -rw-r--r-- | vignettes/web_only/mkin_benchmarks.rda | bin | 854 -> 854 bytes |
10 files changed, 564 insertions, 581 deletions
diff --git a/vignettes/FOCUS_D.Rmd b/vignettes/FOCUS_D.rmd index 02c2c467..b857a3c3 100644 --- a/vignettes/FOCUS_D.Rmd +++ b/vignettes/FOCUS_D.rmd @@ -1,73 +1,73 @@ ----
-title: Example evaluation of FOCUS Example Dataset D
-author: Johannes Ranke
-date: "`r Sys.Date()`"
-output:
-rmarkdown::html_vignette:
- mathjax: null
-vignette: >
- %\VignetteIndexEntry{Example evaluation of FOCUS Example Dataset D}
- %\VignetteEngine{knitr::rmarkdown}
- %\VignetteEncoding{UTF-8}
----
-
-```{r, include = FALSE, message = FALSE}
-library(knitr)
-opts_chunk$set(tidy = FALSE, cache = FALSE)
-library(mkin)
-```
-
-This is just a very simple vignette showing how to fit a degradation model for a parent
-compound with one transformation product using `mkin`. After loading the
-library we look at the data. We have observed concentrations in the column named
-`value` at the times specified in column `time` for the two observed variables
-named `parent` and `m1`.
-
-
-```{r data}
-library(mkin, quietly = TRUE)
-print(FOCUS_2006_D)
-```
-
-Next we specify the degradation model: The parent compound degrades with simple first-order
-kinetics (SFO) to one metabolite named m1, which also degrades with SFO kinetics.
-
-The call to mkinmod returns a degradation model. The differential equations represented in
-R code can be found in the character vector `$diffs` of the `mkinmod` object. If
-a C compiler (gcc) is installed and functional, the differential equation model will
-be compiled from auto-generated C code.
-
-```{r model}
-SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
-print(SFO_SFO$diffs)
-```
-
-We do the fitting without progress report (`quiet = TRUE`).
-
-
-```{r fit}
-fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
-```
-
-A plot of the fit including a residual plot for both observed variables is obtained
-using the `plot_sep` method for `mkinfit` objects, which shows separate graphs for
-all compounds and their residuals.
-
-```{r plot, fig.height = 6, fig.width = 8}
-plot_sep(fit, lpos = c("topright", "bottomright"))
-```
-
-Confidence intervals for the parameter estimates are obtained using the `mkinparplot` function.
-
-
-```{r plot_2, fig.height = 4, fig.width = 8}
-mkinparplot(fit)
-```
-
-A comprehensive report of the results is obtained using the `summary` method for `mkinfit`
-objects.
-
-
-```{r}
-summary(fit)
-```
+--- +title: Example evaluation of FOCUS Example Dataset D +author: Johannes Ranke +date: "`r Sys.Date()`" +output: +rmarkdown::html_vignette: + mathjax: null +vignette: > + %\VignetteIndexEntry{Example evaluation of FOCUS Example Dataset D} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE, message = FALSE} +library(knitr) +opts_chunk$set(tidy = FALSE, cache = FALSE) +library(mkin) +``` + +This is just a very simple vignette showing how to fit a degradation model for a parent +compound with one transformation product using `mkin`. After loading the +library we look at the data. We have observed concentrations in the column named +`value` at the times specified in column `time` for the two observed variables +named `parent` and `m1`. + + +```{r data} +library(mkin, quietly = TRUE) +print(FOCUS_2006_D) +``` + +Next we specify the degradation model: The parent compound degrades with simple first-order +kinetics (SFO) to one metabolite named m1, which also degrades with SFO kinetics. + +The call to mkinmod returns a degradation model. The differential equations represented in +R code can be found in the character vector `$diffs` of the `mkinmod` object. If +a C compiler (gcc) is installed and functional, the differential equation model will +be compiled from auto-generated C code. + +```{r model} +SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +print(SFO_SFO$diffs) +``` + +We do the fitting without progress report (`quiet = TRUE`). + + +```{r fit} +fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) +``` + +A plot of the fit including a residual plot for both observed variables is obtained +using the `plot_sep` method for `mkinfit` objects, which shows separate graphs for +all compounds and their residuals. + +```{r plot, fig.height = 6, fig.width = 8} +plot_sep(fit, lpos = c("topright", "bottomright")) +``` + +Confidence intervals for the parameter estimates are obtained using the `mkinparplot` function. + + +```{r plot_2, fig.height = 4, fig.width = 8} +mkinparplot(fit) +``` + +A comprehensive report of the results is obtained using the `summary` method for `mkinfit` +objects. + + +```{r} +summary(fit) +``` diff --git a/vignettes/FOCUS_L.Rmd b/vignettes/FOCUS_L.rmd index fa6155d2..e017b674 100644 --- a/vignettes/FOCUS_L.Rmd +++ b/vignettes/FOCUS_L.rmd @@ -1,271 +1,271 @@ ----
-title: "Example evaluation of FOCUS Laboratory Data L1 to L3"
-author: "Johannes Ranke"
-date: "`r Sys.Date()`"
-output:
- html_document:
- toc: true
- toc_float:
- collapsed: false
- mathjax: null
- fig_retina: null
-references:
-- id: ranke2014
- title: <span class="nocase">Prüfung und Validierung von Modellierungssoftware als Alternative zu
- ModelMaker 4.0</span>
- author:
- - family: Ranke
- given: Johannes
- type: report
- issued:
- year: 2014
- number: "Umweltbundesamt Projektnummer 27452"
-vignette: >
- %\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3}
- %\VignetteEngine{knitr::rmarkdown}
- %\VignetteEncoding{UTF-8}
----
-
-```{r, include = FALSE}
-library(knitr)
-opts_chunk$set(tidy = FALSE, cache = FALSE)
-```
-
-# Laboratory Data L1
-
-The following code defines example dataset L1 from the FOCUS kinetics
-report, p. 284:
-
-```{r}
-library("mkin", quietly = TRUE)
-FOCUS_2006_L1 = data.frame(
- t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2),
- parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6,
- 72.0, 71.9, 50.3, 59.4, 47.0, 45.1,
- 27.7, 27.3, 10.0, 10.4, 2.9, 4.0))
-FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1)
-```
-
-Here we use the assumptions of simple first order (SFO), the case of declining
-rate constant over time (FOMC) and the case of two different phases of the
-kinetics (DFOP). For a more detailed discussion of the models, please see the
-FOCUS kinetics report.
-
-Since mkin version 0.9-32 (July 2014), we can use shorthand notation like `"SFO"`
-for parent only degradation models. The following two lines fit the model and
-produce the summary report of the model fit. This covers the numerical analysis
-given in the FOCUS report.
-
-```{r}
-m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE)
-summary(m.L1.SFO)
-```
-
-A plot of the fit is obtained with the plot function for mkinfit objects.
-
-```{r fig.width = 6, fig.height = 5}
-plot(m.L1.SFO, show_errmin = TRUE, main = "FOCUS L1 - SFO")
-```
-
-The residual plot can be easily obtained by
-
-```{r fig.width = 6, fig.height = 5}
-mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time")
-```
-
-For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level
-is checked.
-
-```{r fig.width = 6, fig.height = 5}
-m.L1.FOMC <- mkinfit("FOMC", FOCUS_2006_L1_mkin, quiet=TRUE)
-plot(m.L1.FOMC, show_errmin = TRUE, main = "FOCUS L1 - FOMC")
-summary(m.L1.FOMC, data = FALSE)
-```
-
-We get a warning that the default optimisation algorithm `Port` did not converge, which
-is an indication that the model is overparameterised, *i.e.* contains too many
-parameters that are ill-defined as a consequence.
-
-And in fact, due to the higher number of parameters, and the lower number of
-degrees of freedom of the fit, the $\chi^2$ error level is actually higher for
-the FOMC model (3.6%) than for the SFO model (3.4%). Additionally, the
-parameters `log_alpha` and `log_beta` internally fitted in the model have
-excessive confidence intervals, that span more than 25 orders of magnitude (!)
-when backtransformed to the scale of `alpha` and `beta`. Also, the t-test
-for significant difference from zero does not indicate such a significant difference,
-with p-values greater than 0.1, and finally, the parameter correlation of `log_alpha`
-and `log_beta` is 1.000, clearly indicating that the model is overparameterised.
-
-The $\chi^2$ error levels reported in Appendix 3 and Appendix 7 to the FOCUS
-kinetics report are rounded to integer percentages and partly deviate by one
-percentage point from the results calculated by mkin. The reason for
-this is not known. However, mkin gives the same $\chi^2$ error levels
-as the kinfit package and the calculation routines of the kinfit package have
-been extensively compared to the results obtained by the KinGUI
-software, as documented in the kinfit package vignette. KinGUI was the first
-widely used standard package in this field. Also, the calculation of
-$\chi^2$ error levels was compared with KinGUII, CAKE and DegKin manager in
-a project sponsored by the German Umweltbundesamt [@ranke2014].
-
-# Laboratory Data L2
-
-The following code defines example dataset L2 from the FOCUS kinetics
-report, p. 287:
-
-```{r}
-FOCUS_2006_L2 = data.frame(
- t = rep(c(0, 1, 3, 7, 14, 28), each = 2),
- parent = c(96.1, 91.8, 41.4, 38.7,
- 19.3, 22.3, 4.6, 4.6,
- 2.6, 1.2, 0.3, 0.6))
-FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2)
-```
-
-## SFO fit for L2
-
-Again, the SFO model is fitted and the result is plotted. The residual plot
-can be obtained simply by adding the argument `show_residuals` to the plot
-command.
-
-```{r fig.width = 7, fig.height = 6}
-m.L2.SFO <- mkinfit("SFO", FOCUS_2006_L2_mkin, quiet=TRUE)
-plot(m.L2.SFO, show_residuals = TRUE, show_errmin = TRUE,
- main = "FOCUS L2 - SFO")
-```
-
-The $\chi^2$ error level of 14% suggests that the model does not fit very well.
-This is also obvious from the plots of the fit, in which we have included
-the residual plot.
-
-In the FOCUS kinetics report, it is stated that there is no apparent systematic
-error observed from the residual plot up to the measured DT90 (approximately at
-day 5), and there is an underestimation beyond that point.
-
-We may add that it is difficult to judge the random nature of the residuals just
-from the three samplings at days 0, 1 and 3. Also, it is not clear _a
-priori_ why a consistent underestimation after the approximate DT90 should be
-irrelevant. However, this can be rationalised by the fact that the FOCUS fate
-models generally only implement SFO kinetics.
-
-## FOMC fit for L2
-
-For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level
-is checked.
-
-```{r fig.width = 7, fig.height = 6}
-m.L2.FOMC <- mkinfit("FOMC", FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.FOMC, show_residuals = TRUE,
- main = "FOCUS L2 - FOMC")
-summary(m.L2.FOMC, data = FALSE)
-```
-
-The error level at which the $\chi^2$ test passes is much lower in this case.
-Therefore, the FOMC model provides a better description of the data, as less
-experimental error has to be assumed in order to explain the data.
-
-## DFOP fit for L2
-
-Fitting the four parameter DFOP model further reduces the $\chi^2$ error level.
-
-```{r fig.width = 7, fig.height = 6}
-m.L2.DFOP <- mkinfit("DFOP", FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
- main = "FOCUS L2 - DFOP")
-summary(m.L2.DFOP, data = FALSE)
-```
-
-Here, the DFOP model is clearly the best-fit model for dataset L2 based on the
-chi^2 error level criterion. However, the failure to calculate the covariance
-matrix indicates that the parameter estimates correlate excessively. Therefore,
-the FOMC model may be preferred for this dataset.
-
-# Laboratory Data L3
-
-The following code defines example dataset L3 from the FOCUS kinetics report,
-p. 290.
-
-```{r}
-FOCUS_2006_L3 = data.frame(
- t = c(0, 3, 7, 14, 30, 60, 91, 120),
- parent = c(97.8, 60, 51, 43, 35, 22, 15, 12))
-FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3)
-```
-
-## Fit multiple models
-
-As of mkin version 0.9-39 (June 2015), we can fit several models to
-one or more datasets in one call to the function `mmkin`. The datasets
-have to be passed in a list, in this case a named list holding only
-the L3 dataset prepared above.
-
-```{r fig.height = 8}
-# Only use one core here, not to offend the CRAN checks
-mm.L3 <- mmkin(c("SFO", "FOMC", "DFOP"), cores = 1,
- list("FOCUS L3" = FOCUS_2006_L3_mkin), quiet = TRUE)
-plot(mm.L3)
-```
-
-The $\chi^2$ error level of 21% as well as the plot suggest that the SFO model
-does not fit very well. The FOMC model performs better, with an
-error level at which the $\chi^2$ test passes of 7%. Fitting the four
-parameter DFOP model further reduces the $\chi^2$ error level
-considerably.
-
-## Accessing mmkin objects
-
-The objects returned by mmkin are arranged like a matrix, with
-models as a row index and datasets as a column index.
-
-We can extract the summary and plot for *e.g.* the DFOP fit,
-using square brackets for indexing which will result in the use of
-the summary and plot functions working on mkinfit objects.
-
-```{r fig.height = 5}
-summary(mm.L3[["DFOP", 1]])
-plot(mm.L3[["DFOP", 1]], show_errmin = TRUE)
-```
-
-Here, a look to the model plot, the confidence intervals of the parameters
-and the correlation matrix suggest that the parameter estimates are reliable, and
-the DFOP model can be used as the best-fit model based on the $\chi^2$ error
-level criterion for laboratory data L3.
-
-This is also an example where the standard t-test for the parameter `g_ilr` is
-misleading, as it tests for a significant difference from zero. In this case,
-zero appears to be the correct value for this parameter, and the confidence
-interval for the backtransformed parameter `g` is quite narrow.
-
-# Laboratory Data L4
-
-The following code defines example dataset L4 from the FOCUS kinetics
-report, p. 293:
-
-```{r}
-FOCUS_2006_L4 = data.frame(
- t = c(0, 3, 7, 14, 30, 60, 91, 120),
- parent = c(96.6, 96.3, 94.3, 88.8, 74.9, 59.9, 53.5, 49.0))
-FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4)
-```
-
-Fits of the SFO and FOMC models, plots and summaries are produced below:
-
-```{r fig.height = 6}
-# Only use one core here, not to offend the CRAN checks
-mm.L4 <- mmkin(c("SFO", "FOMC"), cores = 1,
- list("FOCUS L4" = FOCUS_2006_L4_mkin),
- quiet = TRUE)
-plot(mm.L4)
-```
-
-The $\chi^2$ error level of 3.3% as well as the plot suggest that the SFO model
-fits very well. The error level at which the $\chi^2$ test passes is slightly
-lower for the FOMC model. However, the difference appears negligible.
-
-
-```{r fig.height = 8}
-summary(mm.L4[["SFO", 1]], data = FALSE)
-summary(mm.L4[["FOMC", 1]], data = FALSE)
-```
-
-
-# References
+--- +title: "Example evaluation of FOCUS Laboratory Data L1 to L3" +author: "Johannes Ranke" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: + collapsed: false + mathjax: null + fig_retina: null +references: +- id: ranke2014 + title: <span class="nocase">Prüfung und Validierung von Modellierungssoftware als Alternative zu + ModelMaker 4.0</span> + author: + - family: Ranke + given: Johannes + type: report + issued: + year: 2014 + number: "Umweltbundesamt Projektnummer 27452" +vignette: > + %\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +library(knitr) +opts_chunk$set(tidy = FALSE, cache = FALSE) +``` + +# Laboratory Data L1 + +The following code defines example dataset L1 from the FOCUS kinetics +report, p. 284: + +```{r} +library("mkin", quietly = TRUE) +FOCUS_2006_L1 = data.frame( + t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2), + parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6, + 72.0, 71.9, 50.3, 59.4, 47.0, 45.1, + 27.7, 27.3, 10.0, 10.4, 2.9, 4.0)) +FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1) +``` + +Here we use the assumptions of simple first order (SFO), the case of declining +rate constant over time (FOMC) and the case of two different phases of the +kinetics (DFOP). For a more detailed discussion of the models, please see the +FOCUS kinetics report. + +Since mkin version 0.9-32 (July 2014), we can use shorthand notation like `"SFO"` +for parent only degradation models. The following two lines fit the model and +produce the summary report of the model fit. This covers the numerical analysis +given in the FOCUS report. + +```{r} +m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE) +summary(m.L1.SFO) +``` + +A plot of the fit is obtained with the plot function for mkinfit objects. + +```{r fig.width = 6, fig.height = 5} +plot(m.L1.SFO, show_errmin = TRUE, main = "FOCUS L1 - SFO") +``` + +The residual plot can be easily obtained by + +```{r fig.width = 6, fig.height = 5} +mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time") +``` + +For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level +is checked. + +```{r fig.width = 6, fig.height = 5} +m.L1.FOMC <- mkinfit("FOMC", FOCUS_2006_L1_mkin, quiet=TRUE) +plot(m.L1.FOMC, show_errmin = TRUE, main = "FOCUS L1 - FOMC") +summary(m.L1.FOMC, data = FALSE) +``` + +We get a warning that the default optimisation algorithm `Port` did not converge, which +is an indication that the model is overparameterised, *i.e.* contains too many +parameters that are ill-defined as a consequence. + +And in fact, due to the higher number of parameters, and the lower number of +degrees of freedom of the fit, the $\chi^2$ error level is actually higher for +the FOMC model (3.6%) than for the SFO model (3.4%). Additionally, the +parameters `log_alpha` and `log_beta` internally fitted in the model have +excessive confidence intervals, that span more than 25 orders of magnitude (!) +when backtransformed to the scale of `alpha` and `beta`. Also, the t-test +for significant difference from zero does not indicate such a significant difference, +with p-values greater than 0.1, and finally, the parameter correlation of `log_alpha` +and `log_beta` is 1.000, clearly indicating that the model is overparameterised. + +The $\chi^2$ error levels reported in Appendix 3 and Appendix 7 to the FOCUS +kinetics report are rounded to integer percentages and partly deviate by one +percentage point from the results calculated by mkin. The reason for +this is not known. However, mkin gives the same $\chi^2$ error levels +as the kinfit package and the calculation routines of the kinfit package have +been extensively compared to the results obtained by the KinGUI +software, as documented in the kinfit package vignette. KinGUI was the first +widely used standard package in this field. Also, the calculation of +$\chi^2$ error levels was compared with KinGUII, CAKE and DegKin manager in +a project sponsored by the German Umweltbundesamt [@ranke2014]. + +# Laboratory Data L2 + +The following code defines example dataset L2 from the FOCUS kinetics +report, p. 287: + +```{r} +FOCUS_2006_L2 = data.frame( + t = rep(c(0, 1, 3, 7, 14, 28), each = 2), + parent = c(96.1, 91.8, 41.4, 38.7, + 19.3, 22.3, 4.6, 4.6, + 2.6, 1.2, 0.3, 0.6)) +FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2) +``` + +## SFO fit for L2 + +Again, the SFO model is fitted and the result is plotted. The residual plot +can be obtained simply by adding the argument `show_residuals` to the plot +command. + +```{r fig.width = 7, fig.height = 6} +m.L2.SFO <- mkinfit("SFO", FOCUS_2006_L2_mkin, quiet=TRUE) +plot(m.L2.SFO, show_residuals = TRUE, show_errmin = TRUE, + main = "FOCUS L2 - SFO") +``` + +The $\chi^2$ error level of 14% suggests that the model does not fit very well. +This is also obvious from the plots of the fit, in which we have included +the residual plot. + +In the FOCUS kinetics report, it is stated that there is no apparent systematic +error observed from the residual plot up to the measured DT90 (approximately at +day 5), and there is an underestimation beyond that point. + +We may add that it is difficult to judge the random nature of the residuals just +from the three samplings at days 0, 1 and 3. Also, it is not clear _a +priori_ why a consistent underestimation after the approximate DT90 should be +irrelevant. However, this can be rationalised by the fact that the FOCUS fate +models generally only implement SFO kinetics. + +## FOMC fit for L2 + +For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level +is checked. + +```{r fig.width = 7, fig.height = 6} +m.L2.FOMC <- mkinfit("FOMC", FOCUS_2006_L2_mkin, quiet = TRUE) +plot(m.L2.FOMC, show_residuals = TRUE, + main = "FOCUS L2 - FOMC") +summary(m.L2.FOMC, data = FALSE) +``` + +The error level at which the $\chi^2$ test passes is much lower in this case. +Therefore, the FOMC model provides a better description of the data, as less +experimental error has to be assumed in order to explain the data. + +## DFOP fit for L2 + +Fitting the four parameter DFOP model further reduces the $\chi^2$ error level. + +```{r fig.width = 7, fig.height = 6} +m.L2.DFOP <- mkinfit("DFOP", FOCUS_2006_L2_mkin, quiet = TRUE) +plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE, + main = "FOCUS L2 - DFOP") +summary(m.L2.DFOP, data = FALSE) +``` + +Here, the DFOP model is clearly the best-fit model for dataset L2 based on the +chi^2 error level criterion. However, the failure to calculate the covariance +matrix indicates that the parameter estimates correlate excessively. Therefore, +the FOMC model may be preferred for this dataset. + +# Laboratory Data L3 + +The following code defines example dataset L3 from the FOCUS kinetics report, +p. 290. + +```{r} +FOCUS_2006_L3 = data.frame( + t = c(0, 3, 7, 14, 30, 60, 91, 120), + parent = c(97.8, 60, 51, 43, 35, 22, 15, 12)) +FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3) +``` + +## Fit multiple models + +As of mkin version 0.9-39 (June 2015), we can fit several models to +one or more datasets in one call to the function `mmkin`. The datasets +have to be passed in a list, in this case a named list holding only +the L3 dataset prepared above. + +```{r fig.height = 8} +# Only use one core here, not to offend the CRAN checks +mm.L3 <- mmkin(c("SFO", "FOMC", "DFOP"), cores = 1, + list("FOCUS L3" = FOCUS_2006_L3_mkin), quiet = TRUE) +plot(mm.L3) +``` + +The $\chi^2$ error level of 21% as well as the plot suggest that the SFO model +does not fit very well. The FOMC model performs better, with an +error level at which the $\chi^2$ test passes of 7%. Fitting the four +parameter DFOP model further reduces the $\chi^2$ error level +considerably. + +## Accessing mmkin objects + +The objects returned by mmkin are arranged like a matrix, with +models as a row index and datasets as a column index. + +We can extract the summary and plot for *e.g.* the DFOP fit, +using square brackets for indexing which will result in the use of +the summary and plot functions working on mkinfit objects. + +```{r fig.height = 5} +summary(mm.L3[["DFOP", 1]]) +plot(mm.L3[["DFOP", 1]], show_errmin = TRUE) +``` + +Here, a look to the model plot, the confidence intervals of the parameters +and the correlation matrix suggest that the parameter estimates are reliable, and +the DFOP model can be used as the best-fit model based on the $\chi^2$ error +level criterion for laboratory data L3. + +This is also an example where the standard t-test for the parameter `g_ilr` is +misleading, as it tests for a significant difference from zero. In this case, +zero appears to be the correct value for this parameter, and the confidence +interval for the backtransformed parameter `g` is quite narrow. + +# Laboratory Data L4 + +The following code defines example dataset L4 from the FOCUS kinetics +report, p. 293: + +```{r} +FOCUS_2006_L4 = data.frame( + t = c(0, 3, 7, 14, 30, 60, 91, 120), + parent = c(96.6, 96.3, 94.3, 88.8, 74.9, 59.9, 53.5, 49.0)) +FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4) +``` + +Fits of the SFO and FOMC models, plots and summaries are produced below: + +```{r fig.height = 6} +# Only use one core here, not to offend the CRAN checks +mm.L4 <- mmkin(c("SFO", "FOMC"), cores = 1, + list("FOCUS L4" = FOCUS_2006_L4_mkin), + quiet = TRUE) +plot(mm.L4) +``` + +The $\chi^2$ error level of 3.3% as well as the plot suggest that the SFO model +fits very well. The error level at which the $\chi^2$ test passes is slightly +lower for the FOMC model. However, the difference appears negligible. + + +```{r fig.height = 8} +summary(mm.L4[["SFO", 1]], data = FALSE) +summary(mm.L4[["FOMC", 1]], data = FALSE) +``` + + +# References diff --git a/vignettes/mkin.Rmd b/vignettes/mkin.rmd index acca0e44..acca0e44 100644 --- a/vignettes/mkin.Rmd +++ b/vignettes/mkin.rmd diff --git a/vignettes/twa.Rmd b/vignettes/twa.rmd index 6f283eb2..6f283eb2 100644 --- a/vignettes/twa.Rmd +++ b/vignettes/twa.rmd diff --git a/vignettes/web_only/FOCUS_Z.Rmd b/vignettes/web_only/FOCUS_Z.rmd index 2da7fde7..2da7fde7 100644 --- a/vignettes/web_only/FOCUS_Z.Rmd +++ b/vignettes/web_only/FOCUS_Z.rmd diff --git a/vignettes/web_only/NAFTA_examples.Rmd b/vignettes/web_only/NAFTA_examples.rmd index 26a9240a..26a9240a 100644 --- a/vignettes/web_only/NAFTA_examples.Rmd +++ b/vignettes/web_only/NAFTA_examples.rmd 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 @@ <meta name="author" content="Johannes Ranke" /> -<meta name="date" content="2020-05-12" /> +<meta name="date" content="2020-05-13" /> <title>Benchmark timings for mkin</title> @@ -1583,29 +1583,12 @@ div.tocify { <h1 class="title toc-ignore">Benchmark timings for mkin</h1> <h4 class="author">Johannes Ranke</h4> -<h4 class="date">2020-05-12</h4> +<h4 class="date">2020-05-13</h4> </div> -<p>Each system is characterized by its CPU type, the operating system type and the mkin version. Currently only values for one system are available.</p> -<pre class="r"><code>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) - } -}</code></pre> +<p>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.</p> <div id="test-cases" class="section level2"> <h2>Test cases</h2> <p>Parent only:</p> @@ -1619,17 +1602,14 @@ t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP <p>One metabolite:</p> <pre class="r"><code>SFO_SFO <- mkinmod( parent = mkinsub("SFO", "m1"), - m1 = mkinsub("SFO"))</code></pre> -<pre><code>## Successfully compiled differential equation model from auto-generated C code.</code></pre> -<pre class="r"><code>FOMC_SFO <- mkinmod( + m1 = mkinsub("SFO")) +FOMC_SFO <- mkinmod( parent = mkinsub("FOMC", "m1"), - m1 = mkinsub("SFO"))</code></pre> -<pre><code>## Successfully compiled differential equation model from auto-generated C code.</code></pre> -<pre class="r"><code>DFOP_SFO <- mkinmod( + m1 = mkinsub("SFO")) +DFOP_SFO <- mkinmod( parent = mkinsub("FOMC", "m1"), - m1 = mkinsub("SFO"))</code></pre> -<pre><code>## Successfully compiled differential equation model from auto-generated C code.</code></pre> -<pre class="r"><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 </div> <div id="results" class="section level2"> <h2>Results</h2> -<p>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.</p> -<pre class="r"><code>rownames(mkin_benchmarks) <- as.character(mkin_benchmarks$mkin)</code></pre> +<p>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.</p> <p>Benchmarks for all available error models are shown.</p> <div id="parent-only" class="section level3"> <h3>Parent only</h3> -<p>Constant variance and two-component error model:</p> -<pre class="r"><code>kable(mkin_benchmarks[, c("t1", "t2")])</code></pre> +<p>Constant variance (t1) and two-component error model (t2) for four models fitted to two datasets, i.e. eight fits for each test.</p> <table> <thead> <tr class="header"> -<th></th> -<th align="right">t1</th> -<th align="right">t2</th> +<th align="left">mkin version</th> +<th align="right">t1 [s]</th> +<th align="right">t2 [s]</th> </tr> </thead> <tbody> <tr class="odd"> -<td>0.9.48.1</td> +<td align="left">0.9.48.1</td> <td align="right">3.610</td> <td align="right">11.019</td> </tr> <tr class="even"> -<td>0.9.49.1</td> +<td align="left">0.9.49.1</td> <td align="right">8.184</td> <td align="right">22.889</td> </tr> <tr class="odd"> -<td>0.9.49.2</td> +<td align="left">0.9.49.2</td> <td align="right">7.064</td> <td align="right">12.558</td> </tr> <tr class="even"> -<td>0.9.49.3</td> +<td align="left">0.9.49.3</td> <td align="right">7.296</td> <td align="right">21.239</td> </tr> <tr class="odd"> -<td>0.9.49.4</td> +<td align="left">0.9.49.4</td> <td align="right">5.936</td> <td align="right">20.545</td> </tr> <tr class="even"> -<td>0.9.50.2</td> -<td align="right">1.559</td> -<td align="right">3.929</td> +<td align="left">0.9.50.2</td> +<td align="right">1.547</td> +<td align="right">3.928</td> </tr> </tbody> </table> </div> <div id="one-metabolite" class="section level3"> <h3>One metabolite</h3> -<p>Constant variance, variance by variable and two-component error model:</p> -<pre class="r"><code>kable(mkin_benchmarks[, c("t3", "t4", "t5")])</code></pre> +<p>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.</p> <table> <thead> <tr class="header"> -<th></th> -<th align="right">t3</th> -<th align="right">t4</th> -<th align="right">t5</th> +<th align="left">mkin version</th> +<th align="right">t3 [s]</th> +<th align="right">t4 [s]</th> +<th align="right">t5 [s]</th> </tr> </thead> <tbody> <tr class="odd"> -<td>0.9.48.1</td> +<td align="left">0.9.48.1</td> <td align="right">3.764</td> <td align="right">14.347</td> <td align="right">9.495</td> </tr> <tr class="even"> -<td>0.9.49.1</td> +<td align="left">0.9.49.1</td> <td align="right">4.649</td> <td align="right">13.789</td> <td align="right">6.395</td> </tr> <tr class="odd"> -<td>0.9.49.2</td> +<td align="left">0.9.49.2</td> <td align="right">4.786</td> <td align="right">8.461</td> <td align="right">5.675</td> </tr> <tr class="even"> -<td>0.9.49.3</td> +<td align="left">0.9.49.3</td> <td align="right">4.510</td> <td align="right">13.805</td> <td align="right">7.386</td> </tr> <tr class="odd"> -<td>0.9.49.4</td> +<td align="left">0.9.49.4</td> <td align="right">4.446</td> <td align="right">15.335</td> <td align="right">6.002</td> </tr> <tr class="even"> -<td>0.9.50.2</td> -<td align="right">1.352</td> -<td align="right">6.110</td> -<td align="right">2.841</td> +<td align="left">0.9.50.2</td> +<td align="right">1.371</td> +<td align="right">6.154</td> +<td align="right">2.720</td> </tr> </tbody> </table> </div> <div id="two-metabolites" class="section level3"> <h3>Two metabolites</h3> -<p>Two different datasets, for each constant variance, variance by variable and two-component error model are shown:</p> -<pre class="r"><code>kable(mkin_benchmarks[, paste0("t", 6:11)])</code></pre> +<p>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.</p> <table> <thead> <tr class="header"> -<th></th> -<th align="right">t6</th> -<th align="right">t7</th> -<th align="right">t8</th> -<th align="right">t9</th> -<th align="right">t10</th> -<th align="right">t11</th> +<th align="left">mkin version</th> +<th align="right">t6 [s]</th> +<th align="right">t7 [s]</th> +<th align="right">t8 [s]</th> +<th align="right">t9 [s]</th> +<th align="right">t10 [s]</th> +<th align="right">t11 [s]</th> </tr> </thead> <tbody> <tr class="odd"> -<td>0.9.48.1</td> +<td align="left">0.9.48.1</td> <td align="right">2.623</td> <td align="right">4.587</td> <td align="right">7.525</td> @@ -1796,7 +1772,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks <td align="right">31.267</td> </tr> <tr class="even"> -<td>0.9.49.1</td> +<td align="left">0.9.49.1</td> <td align="right">2.542</td> <td align="right">4.128</td> <td align="right">4.632</td> @@ -1805,7 +1781,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks <td align="right">5.636</td> </tr> <tr class="odd"> -<td>0.9.49.2</td> +<td align="left">0.9.49.2</td> <td align="right">2.723</td> <td align="right">4.478</td> <td align="right">4.862</td> @@ -1814,7 +1790,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks <td align="right">5.574</td> </tr> <tr class="even"> -<td>0.9.49.3</td> +<td align="left">0.9.49.3</td> <td align="right">2.643</td> <td align="right">4.374</td> <td align="right">7.020</td> @@ -1823,7 +1799,7 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks <td align="right">7.365</td> </tr> <tr class="odd"> -<td>0.9.49.4</td> +<td align="left">0.9.49.4</td> <td align="right">2.635</td> <td align="right">4.259</td> <td align="right">4.737</td> @@ -1832,13 +1808,13 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks <td align="right">5.626</td> </tr> <tr class="even"> -<td>0.9.50.2</td> -<td align="right">0.759</td> -<td align="right">1.204</td> -<td align="right">1.275</td> -<td align="right">2.837</td> -<td align="right">2.026</td> -<td align="right">2.976</td> +<td align="left">0.9.50.2</td> +<td align="right">0.742</td> +<td align="right">1.203</td> +<td align="right">1.282</td> +<td align="right">2.845</td> +<td align="right">2.037</td> +<td align="right">2.987</td> </tr> </tbody> </table> diff --git a/vignettes/web_only/benchmarks.Rmd b/vignettes/web_only/benchmarks.rmd index 990c2fee..7ae12451 100644 --- a/vignettes/web_only/benchmarks.Rmd +++ b/vignettes/web_only/benchmarks.rmd @@ -20,9 +20,11 @@ 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. +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} +```{r include = FALSE} cpu_model <- benchmarkme::get_cpu()$model_name operating_system <- Sys.info()[["sysname"]] mkin_version <- as.character(packageVersion("mkin")) @@ -58,7 +60,7 @@ t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), parent_datasets, One metabolite: -```{r one_metabolite} +```{r one_metabolite, message = FALSE} SFO_SFO <- mkinmod( parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) @@ -77,7 +79,7 @@ t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_D), Two metabolites, synthetic data: -```{r two_metabolites} +```{r two_metabolites, message = FALSE} m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"), M1 = mkinsub("SFO", "M2"), M2 = mkinsub("SFO"), @@ -116,9 +118,11 @@ save(mkin_benchmarks, file = "~/git/mkin/vignettes/web_only/mkin_benchmarks.rda" 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} -rownames(mkin_benchmarks) <- as.character(mkin_benchmarks$mkin) +```{r, include = FALSE} +dimnames(mkin_benchmarks) <- list(as.character(mkin_benchmarks$mkin), + c("CPU", "OS", "mkin version", paste0("t", 1:11, " [s]"))) ``` @@ -126,25 +130,28 @@ Benchmarks for all available error models are shown. ### Parent only -Constant variance and two-component error model: +Constant variance (t1) and two-component error model (t2) for four models +fitted to two datasets, i.e. eight fits for each test. -```{r} -kable(mkin_benchmarks[, c("t1", "t2")]) +```{r, echo = FALSE} +kable(mkin_benchmarks[, 4:5], rownames = "mkin version") ``` ### One metabolite -Constant variance, variance by variable and two-component error model: +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} -kable(mkin_benchmarks[, c("t3", "t4", "t5")]) +```{r, echo = FALSE} +kable(mkin_benchmarks[, 6:8], rownames = "mkin version") ``` ### Two metabolites -Two different datasets, for each constant variance, variance by variable and -two-component error model are shown: +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} -kable(mkin_benchmarks[, paste0("t", 6:11)]) +```{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 index f99ea808..0b8c617a 100644 --- a/vignettes/web_only/compiled_models.Rmd +++ b/vignettes/web_only/compiled_models.rmd @@ -1,141 +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]]))
-}
-```
+--- +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 Binary files differindex 268e5efe..c88ca4d6 100644 --- a/vignettes/web_only/mkin_benchmarks.rda +++ b/vignettes/web_only/mkin_benchmarks.rda |