diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | DESCRIPTION | 14 | ||||
| -rw-r--r-- | R/hierarchical_kinetics.R | 1 | ||||
| -rw-r--r-- | inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd | 181 |
4 files changed, 175 insertions, 22 deletions
@@ -3,6 +3,7 @@ install.log inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_cache inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_files inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.pdf +inst/rmarkdown/templates/hierarchical_kinetics/skeleton/*_dlls mkin*.tar.gz mkin.tar mkin.Rcheck diff --git a/DESCRIPTION b/DESCRIPTION index 04237b10..d061a5a9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data Version: 1.2.2 -Date: 2022-11-22 +Date: 2023-01-03 Authors@R: c( person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "johannes.ranke@jrwb.de", @@ -17,11 +17,11 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, equation models are solved using automatically generated C functions. Heteroscedasticity can be taken into account using variance by variable or two-component error models as described by Ranke and Meinecke (2018) - <doi:10.3390/environments6120124>. Interfaces to several nonlinear - mixed-effects model packages are available, some of which are described by - Ranke et al. (2021) <doi:10.3390/environments8080071>. Please note that no - warranty is implied for correctness of results or fitness for a particular - purpose. + <doi:10.3390/environments6120124>. Hierarchical degradation models can + be fitted using nonlinear mixed-effects model packages as a backend as + described by Ranke et al. (2021) <doi:10.3390/environments8080071>. Please + note that no warranty is implied for correctness of results or fitness for a + particular purpose. Depends: R (>= 2.15.1), Imports: stats, graphics, methods, parallel, deSolve, R6, inline (>= 0.3.19), numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), saemix (>= 3.1), rlang, vctrs @@ -36,4 +36,4 @@ VignetteBuilder: knitr BugReports: https://github.com/jranke/mkin/issues/ URL: https://pkgdown.jrwb.de/mkin/ Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.2 +RoxygenNote: 7.2.3 diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R index 70e0100a..a90dd619 100644 --- a/R/hierarchical_kinetics.R +++ b/R/hierarchical_kinetics.R @@ -31,6 +31,7 @@ hierarchical_kinetics <- function(..., keep_tex = FALSE) { fmt <- rmarkdown::pdf_document(..., keep_tex = keep_tex, toc = TRUE, + toc_depth = 3, includes = rmarkdown::includes(in_header = "header.tex"), extra_dependencies = c("float", "listing", "framed") ) diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd index 3c7dc8a5..9a975607 100644 --- a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd @@ -6,6 +6,10 @@ output: mkin::hierarchical_kinetics geometry: margin=2cm --- +\clearpage + +# Setup + ```{r packages, cache = FALSE, message = FALSE} library(mkin) library(knitr) @@ -14,7 +18,7 @@ library(parallel) library(readxl) ``` -```{r n_cores, cache = FALSE, echo = FALSE} +```{r n_cores, cache = FALSE} n_cores <- detectCores() if (Sys.info()["sysname"] == "Windows") { @@ -92,7 +96,8 @@ illparms(parent_mhmkin) |> kable() ``` Therefore, the fits are updated, excluding random effects that were -ill-defined according to the `illparms` function. +ill-defined according to the `illparms` function. The status of the fits +is checked. ```{r parent-mhmkin-refined} parent_mhmkin_refined <- update(parent_mhmkin, @@ -100,14 +105,22 @@ parent_mhmkin_refined <- update(parent_mhmkin, status(parent_mhmkin_refined) |> kable() ``` -The most suitable model is selected based on the AIC. +Also, it is checked if the AIC values of the refined fits are actually smaller +than the AIC values of the original fits. + +```{r dependson = "parent-mhmkin-refined"} +(AIC(parent_mhmkin_refined) < AIC(parent_mhmkin)) |> kable() +``` + +From the refined fits, the most suitable model is selected using the AIC. -```{r, dependson = "parent-mhmkin"} +```{r parent-best, dependson = "parent-mhmkin"} aic_parent <- AIC(parent_mhmkin_refined) min_aic <- which(aic_parent == min(aic_parent), arr.ind = TRUE) best_degmod_parent <- rownames(aic_parent)[min_aic[1]] best_errmod_parent <- colnames(aic_parent)[min_aic[2]] anova(parent_mhmkin_refined) |> kable(digits = 1) +parent_best <- parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]] ``` Based on the AIC, the combination of the `r best_degmod_parent` degradation @@ -115,25 +128,152 @@ model with the error model `r errmods[best_errmod_parent]` is identified to be most suitable for the degradation of the parent. The check below confirms that no ill-defined parameters remain for this combined model. -```{r dependson = "parent-mhmkin"} -illparms(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]]) +```{r dependson = "parent-best"} +illparms(parent_best) ``` -The corresponding fit is shown below. +The corresponding fit is plotted below. + +```{r dependson = "parent-best"} +plot(parent_best) +``` +The fitted parameters, together with approximate confidence +intervals are listed below. -```{r parent-best-full, dependson = "parent-mhmkin"} -plot(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]]) +```{r dependson = "parent-best"} +parms(parent_best, ci = TRUE) |> kable(digits = 3) ``` -Detailed listings of the parent fits are shown in the Appendix. +As an example of the investigation of covariate influence on degradation +parameters, a covariate model is added to the hierarchical model for each of +the degradation parameters with well-defined random effects. Also, +a version with covariate models for both of them is fitted. + +```{r parent-best-pH} +parent_best_pH_1 <- update(parent_best, covariates = covariates, + covariate_models = list(log_k_lambda_free ~ pH)) +parent_best_pH_2 <- update(parent_best, covariates = covariates, + covariate_models = list(log_k_lambda_bound_free ~ pH)) +parent_best_pH_3 <- update(parent_best, covariates = covariates, + covariate_models = list(log_k_lambda_free ~ pH, log_k_lambda_bound_free ~ pH)) +``` + +The resulting models are compared. + +```{r dependson = "parent-best-pH"} +anova(parent_best, parent_best_pH_1, parent_best_pH_2, parent_best_pH_3) |> + kable(digits = 1) +``` + +The model fit with the lowest AIC is the one with a pH correlation of the +desorption rate constant `k_lambda_bound_free`. Plot and parameter listing +of this fit are shown below. Also, it is confirmed that no ill-defined +variance parameters are found. + +```{r dependson = "parent-best-pH"} +plot(parent_best_pH_2) +``` + +```{r dependson = "parent-best-pH"} +illparms(parent_best_pH_2) +parms(parent_best_pH_2, ci = TRUE) |> kable(digits = 3) +``` \clearpage -# Appendix +# Pathway fits + +As an example of a pathway fit, a model with SFORB for the parent compound and +parallel formation of two metabolites is set up. + + +```{r m_sforb_sfo2} +if (!dir.exists("lambda_dlls")) dir.create("lambda_dlls") + +m_sforb_sfo2 = mkinmod( + lambda = mkinsub("SFORB", to = c("c_V", "c_XV")), + c_V = mkinsub("SFO"), + c_XV = mkinsub("SFO"), + name = "sforb_sfo2", + dll_dir = "lambda_dlls", + overwrite = TRUE, quiet = TRUE +) +``` + +Separate evaluations of all datasets are performed with constant variance +and using two-component error. -## Summaries of saem fits +```{r path-sep, dependson = c("m_sforb_all", "ds")} +sforb_sep_const <- mmkin(list(sforb_path = m_sforb_sfo2), ds, + cluster = cl, quiet = TRUE) +sforb_sep_tc <- update(sforb_sep_const, error_model = "tc") +``` + +The separate fits with constant variance are plotted. + +```{r dependson = "path-sep", fig.height = 9} +plot(mixed(sforb_sep_const)) +``` -### Parent fits +The two corresponding hierarchical fits, with the random effects for the parent +degradation parameters excluded as discussed above, and including the covariate +model that was identified for the parent degradation, are attempted below. + +```{r path-1, dependson = "path-sep"} +path_1 <- mhmkin(list(sforb_sep_const, sforb_sep_tc), + no_random_effect = c("lambda_free_0", "log_k_lambda_free_bound"), + covariates = covariates, covariate_models = list(log_k_lambda_bound_free ~ pH), + cluster = cl) +``` + +```{r dependson = "path-1"} +status(path_1) |> kable() +``` + +The status information shows that both fits were successfully completed. + +```{r dependson = "path-1"} +anova(path_1) |> kable(digits = 1) +``` +Model comparison shows that the two-component error model provides a much +better fit. + +```{r dependson = "path-1"} +illparms(path_1[["sforb_path", "tc"]]) +``` + +Two ill-defined variance components are found. Therefore, the fit is +repeated with the corresponding random effects removed. + +```{r path-1-refined, dependson = "path-1"} +path_1_refined <- update(path_1[["sforb_path", "tc"]], + no_random_effect = c("lambda_free_0", "log_k_lambda_free_bound", + "log_k_c_XV", "f_lambda_ilr_2")) +``` + +The empty output of the illparms function indicates that there are no +ill-defined parameters remaining in the refined fit. + +```{r dependson = "path-1-refined"} +illparms(path_1_refined) +``` + +Below, the refined fit is plotted and the fitted parameters are shown together +with their 95% confidence intervals. + +```{r dependson = "path-1-refined", fig.height = 9} +plot(path_1_refined) +``` + +```{r dependson = "path-1-refined", fig.height = 9} +parms(path_1_refined, ci = TRUE) |> kable(digits = 3) +``` + +\clearpage + +# Appendix + +## Listings of initial parent fits ```{r listings-parent, results = "asis", echo = FALSE, dependson = "parent_mhmkin"} for (deg_mod in parent_deg_mods) { @@ -144,7 +284,7 @@ for (deg_mod in parent_deg_mods) { } ``` -### Refined parent fits +## Listings of refined parent fits ```{r listings-parent-refined, results = "asis", echo = FALSE, dependson = "parent_mhmkin_refined"} for (deg_mod in parent_deg_mods) { @@ -155,9 +295,20 @@ for (deg_mod in parent_deg_mods) { } ``` +## Listings of pathway fits + +```{r listings-path-1, results = "asis", echo = FALSE, dependson = "path-1-refined"} +tex_listing(path_1[["sforb_path", "const"]], + caption = "Hierarchical fit of SFORB-SFO2 with constant variance") +tex_listing(path_1[["sforb_path", "tc"]], + caption = "Hierarchical fit of SFORB-SFO2 with two-component error") +tex_listing(path_1_refined, + caption = "Refined hierarchical fit of SFORB-SFO2 with two-component error") +``` + ## Session info -```{r, echo = FALSE, cache = FALSE} +```{r echo = FALSE, cache = FALSE} parallel::stopCluster(cl) sessionInfo() ``` |