From 0023df3c31fac29b5f9337ecd732a5dfd4d51a2d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 19 Dec 2022 06:37:32 +0100 Subject: Template and spreadsheet for hierarchical kinetics The template only shows parent data evaluation without covariate models for now. The spreadsheet will also be useful for unit testing of the read_spreadsheet function. --- .../rmarkdown/templates/hier/skeleton/skeleton.Rmd | 186 +++++++++++++++++++++ inst/rmarkdown/templates/hier/template.yaml | 3 + 2 files changed, 189 insertions(+) create mode 100644 inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd create mode 100644 inst/rmarkdown/templates/hier/template.yaml (limited to 'inst/rmarkdown') diff --git a/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd new file mode 100644 index 00000000..df6fa2f9 --- /dev/null +++ b/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd @@ -0,0 +1,186 @@ +--- +title: "Hierarchical kinetic modelling of degradation data" +author: +date: Last change on DD MMM YYYY, last compiled on `r format(Sys.time(), + "%e %B %Y")` +output: + pdf_document: + extra_dependencies: ["float", "listing"] +toc: yes +geometry: margin=2cm +--- + + +```{r setup, echo = FALSE, cache = FALSE} +errmods <- c(const = "constant variance", tc = "two-component error") + +knitr::opts_chunk$set( + comment = "", tidy = FALSE, cache = TRUE, fig.pos = "H", fig.align = "center" +) +options(knitr.kable.NA = "") + +# Version requirements +if (getRversion() < "4.1.0") + stop("You need R with version > 4.1.0 to compile this document") +if ((saemix_version <- packageVersion("saemix")) < "3.1") { + warning("Your saemix version is ", saemix_version, + ", you should preferably use 3.2 to compile this document") +} +if ((mkin_version <- packageVersion("mkin")) < "1.2.2") { + stop("Your mkin version is ", mkin_version, + ", you need at least 1.2.2 to compile this document") +} +``` + +```{r packages, cache = FALSE, message = FALSE, warning = FALSE, echo = FALSE} +library(mkin) +library(saemix) +library(parallel) +library(knitr) +``` + +```{r n_cores, cache = FALSE, echo = FALSE} +n_cores <- detectCores() + +if (Sys.info()["sysname"] == "Windows") { + cl <- makePSOCKcluster(n_cores) +} else { + cl <- makeForkCluster(n_cores) +} +``` + +\clearpage + +# Introduction + +This report shows hierarchical kinetic modelling for ... +The data were obtained from ... + +```{r ds} +data_path <- system.file("testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", package = "mkin") +ds <- read_spreadsheet(data_path, valid_datasets = c(1:4, 7:13)) +covariates <- attr(ds, "covariates") +``` + +The covariate data are shown below. + +```{r results = "asis", dependson = "ds", echo = FALSE} +kable(covariates, caption = "Covariate data for all datasets") +``` + +\clearpage + +The datasets with the residue time series are shown in the tables below. Please +refer to the spreadsheet for details like data sources, treatment of values +below reporting limits and time step normalisation factors. + +```{r results = "asis", dependson = "ds", echo = FALSE} +for (ds_name in names(ds)) { + print( + kable(mkin_long_to_wide(ds[[ds_name]]), + caption = paste("Dataset", ds_name), + booktabs = TRUE, row.names = FALSE)) + cat("\n\\clearpage\n") +} +``` + +# Parent only evaluations + +The following code performs separate fits of the candidate degradation models +to all datasets using constant variance and the two-component error model. + +```{r parent-sep, dependson = "ds"} +parent_deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB") +parent_sep_const <- mmkin( + parent_deg_mods, ds, + error_model = "const", + cluster = cl, quiet = TRUE) +parent_sep_tc <- update(parent_sep_const, error_model = "tc") +``` + +To select the parent model, the corresponding hierarchical fits are performed below. + +```{r parent-mhmkin, dependson = "parent-sep"} +parent_mhmkin <- mhmkin(list(parent_sep_const, parent_sep_tc), cluster = cl) +status(parent_mhmkin) |> kable() +``` + +All fits terminate without errors (status OK). The check for ill-defined +parameters shows that not all random effect parameters can be robustly +quantified. + +```{r dependson = "parent_mhmkin"} +illparms(parent_mhmkin) |> kable() +``` + +Therefore, the fits are updated, excluding random effects that were +ill-defined according to the `illparms` function. + +```{r parent-mhmkin-refined} +parent_mhmkin_refined <- update(parent_mhmkin, + no_random_effect = illparms(parent_mhmkin)) +status(parent_mhmkin_refined) |> kable() +``` + +The most suitable model is selected based on the AIC. + +```{r 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) +``` + +Based on the AIC, the combination of the `r best_degmod_parent` degradation +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]]) +``` + +The corresponding fit is shown below. + +```{r parent-best-full, dependson = "parent-mhmkin"} +plot(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]]) +``` + +Detailed listings of the parent fits are shown in the Appendix. + +\clearpage + +# Appendix + +## Summaries of saem fits + +### Parent fits + +```{r listings-parent, results = "asis", echo = FALSE} +for (deg_mod in parent_deg_mods) { + for (err_mod in c("const", "tc")) { + caption <- paste("Hierarchical", deg_mod, "fit with", errmods[err_mod]) + tex_listing(parent_mhmkin[[deg_mod, err_mod]], caption) + } +} +``` + +### Refined parent fits + +```{r listings-pathway, results = "asis", echo = FALSE} +for (deg_mod in parent_deg_mods) { + for (err_mod in c("const", "tc")) { + caption <- paste("Refined hierarchical", deg_mod, "fit with", errmods[err_mod]) + tex_listing(parent_mhmkin_refined[[deg_mod, err_mod]], caption) + } +} +``` + +## Session info + +```{r, echo = FALSE} +parallel::stopCluster(cl) +sessionInfo() +``` + diff --git a/inst/rmarkdown/templates/hier/template.yaml b/inst/rmarkdown/templates/hier/template.yaml new file mode 100644 index 00000000..d8ab6a4d --- /dev/null +++ b/inst/rmarkdown/templates/hier/template.yaml @@ -0,0 +1,3 @@ +name: Hierarchical kinetics +description: Hierarchical kinetic modelling of degradation data +create_dir: true -- cgit v1.2.3 From 886c9ef013124aa954d960c655b349b5340ff154 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 19 Dec 2022 12:31:56 +0100 Subject: Rename template folder, create format Instead of rmarkdown::pdf_document, mkin::hierarchical_kinetics is used as a document format in the template. In this way, the template file can be freed from some R code and yaml options that the average user does not have to be aware of. --- .Rbuildignore | 6 +- .gitignore | 6 +- GNUmakefile | 4 +- NAMESPACE | 2 + NEWS.md | 4 + R/hierarchical_kinetics.R | 39 +++++ R/parplot.R | 2 +- .../rmarkdown/templates/hier/skeleton/skeleton.Rmd | 186 --------------------- inst/rmarkdown/templates/hier/template.yaml | 3 - .../hierarchical_kinetics/skeleton/header.tex | 1 + .../hierarchical_kinetics/skeleton/skeleton.Rmd | 162 ++++++++++++++++++ .../templates/hierarchical_kinetics/template.yaml | 3 + log/check.log | 21 ++- man/hierarchical_kinetics.Rd | 29 ++++ 14 files changed, 265 insertions(+), 203 deletions(-) create mode 100644 R/hierarchical_kinetics.R delete mode 100644 inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd delete mode 100644 inst/rmarkdown/templates/hier/template.yaml create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/template.yaml create mode 100644 man/hierarchical_kinetics.Rd (limited to 'inst/rmarkdown') diff --git a/.Rbuildignore b/.Rbuildignore index 0dcec29a..1eb33577 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,9 +9,9 @@ ^test.R$ ^mkin_.*\.tar\.gz ^mkin.tar$ -^inst/rmarkdown/templates/hier/skeleton/skeleton.pdf$ -^inst/rmarkdown/templates/hier/skeleton/skeleton_cache -^inst/rmarkdown/templates/hier/skeleton/skeleton_files +^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.pdf$ +^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_cache +^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_files ^vignettes/.build.timestamp$ ^vignettes/.*_cache$ ^vignettes/.*cache$ diff --git a/.gitignore b/.gitignore index 4f479259..9808896d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,8 @@ docs/articles/*_cache/ install.log -inst/rmarkdown/templates/hier/skeleton/skeleton_cache -inst/rmarkdown/templates/hier/skeleton/skeleton_files -inst/rmarkdown/templates/hier/skeleton/skeleton.pdf +inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_cache +inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_files +inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.pdf mkin*.tar.gz mkin.tar mkin.Rcheck diff --git a/GNUmakefile b/GNUmakefile index 4680514c..4a33c538 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -24,8 +24,8 @@ pkgfiles = \ DESCRIPTION \ inst/WORDLIST \ inst/dataset_generation/* \ - inst/rmarkdown/templates/hier/template.yaml \ - inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd \ + inst/rmarkdown/templates/hierarchical_kinetics/template.yaml \ + inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd \ inst/testdata/* \ man/* \ NAMESPACE \ diff --git a/NAMESPACE b/NAMESPACE index 107ffc54..84d6b713 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -96,6 +96,7 @@ export(create_deg_func) export(endpoints) export(f_time_norm_focus) export(get_deg_func) +export(hierarchical_kinetics) export(illparms) export(ilr) export(intervals) @@ -187,6 +188,7 @@ importFrom(stats,qf) importFrom(stats,qlogis) importFrom(stats,qnorm) importFrom(stats,qt) +importFrom(stats,quantile) importFrom(stats,residuals) importFrom(stats,rnorm) importFrom(stats,shapiro.test) diff --git a/NEWS.md b/NEWS.md index 7e65204f..2c09df35 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # mkin 1.2.2 +- 'inst/rmarkdown/templates/hier': R markdown template to facilitate the application of hierarchical kinetic models. + +- 'inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx': Example spreadsheet for use with 'read_spreadsheet()'. + - 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well. - 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted. List correlations of random effects if specified by the user in the covariance model. diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R new file mode 100644 index 00000000..f7ffb333 --- /dev/null +++ b/R/hierarchical_kinetics.R @@ -0,0 +1,39 @@ +#' Hierarchical kinetics template +#' +#' R markdown format for setting up hierarchical kinetics based on a template +#' provided with the mkin package. +#' +#' @inheritParams rmarkdown::pdf_document +#' @param ... Arguments to \code{rmarkdown::pdf_document} +#' +#' @return R Markdown output format to pass to +#' \code{\link[rmarkdown:render]{render}} +#' +#' @examples +#' +#' \dontrun{ +#' library(rmarkdown) +#' draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin") +#' } +#' +#' @export +hierarchical_kinetics <- function(..., keep_tex = FALSE) { + + if (getRversion() < "4.1.0") + stop("You need R with version > 4.1.0 to compile this document") + + if (!requireNamespace("knitr")) stop("Please install the knitr package to use this template") + if (!requireNamespace("rmarkdown")) stop("Please install the rmarkdown package to use this template") + knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE) + knitr::opts_chunk$set(fig.align = "center", fig.pos = "H") + options(knitr.kable.NA = "") + + fmt <- rmarkdown::pdf_document(..., + keep_tex = keep_tex, + toc = TRUE, + includes = rmarkdown::includes(in_header = "header.tex"), + extra_dependencies = c("float", "listing", "framed") + ) + + return(fmt) +} diff --git a/R/parplot.R b/R/parplot.R index e9c18947..3da4b51a 100644 --- a/R/parplot.R +++ b/R/parplot.R @@ -23,7 +23,7 @@ #' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. #' doi: 10.1186/s12859-021-04373-4. #' @seealso [multistart] -#' @importFrom stats median +#' @importFrom stats median quantile #' @export parplot <- function(object, ...) { UseMethod("parplot") diff --git a/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd deleted file mode 100644 index df6fa2f9..00000000 --- a/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd +++ /dev/null @@ -1,186 +0,0 @@ ---- -title: "Hierarchical kinetic modelling of degradation data" -author: -date: Last change on DD MMM YYYY, last compiled on `r format(Sys.time(), - "%e %B %Y")` -output: - pdf_document: - extra_dependencies: ["float", "listing"] -toc: yes -geometry: margin=2cm ---- - - -```{r setup, echo = FALSE, cache = FALSE} -errmods <- c(const = "constant variance", tc = "two-component error") - -knitr::opts_chunk$set( - comment = "", tidy = FALSE, cache = TRUE, fig.pos = "H", fig.align = "center" -) -options(knitr.kable.NA = "") - -# Version requirements -if (getRversion() < "4.1.0") - stop("You need R with version > 4.1.0 to compile this document") -if ((saemix_version <- packageVersion("saemix")) < "3.1") { - warning("Your saemix version is ", saemix_version, - ", you should preferably use 3.2 to compile this document") -} -if ((mkin_version <- packageVersion("mkin")) < "1.2.2") { - stop("Your mkin version is ", mkin_version, - ", you need at least 1.2.2 to compile this document") -} -``` - -```{r packages, cache = FALSE, message = FALSE, warning = FALSE, echo = FALSE} -library(mkin) -library(saemix) -library(parallel) -library(knitr) -``` - -```{r n_cores, cache = FALSE, echo = FALSE} -n_cores <- detectCores() - -if (Sys.info()["sysname"] == "Windows") { - cl <- makePSOCKcluster(n_cores) -} else { - cl <- makeForkCluster(n_cores) -} -``` - -\clearpage - -# Introduction - -This report shows hierarchical kinetic modelling for ... -The data were obtained from ... - -```{r ds} -data_path <- system.file("testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", package = "mkin") -ds <- read_spreadsheet(data_path, valid_datasets = c(1:4, 7:13)) -covariates <- attr(ds, "covariates") -``` - -The covariate data are shown below. - -```{r results = "asis", dependson = "ds", echo = FALSE} -kable(covariates, caption = "Covariate data for all datasets") -``` - -\clearpage - -The datasets with the residue time series are shown in the tables below. Please -refer to the spreadsheet for details like data sources, treatment of values -below reporting limits and time step normalisation factors. - -```{r results = "asis", dependson = "ds", echo = FALSE} -for (ds_name in names(ds)) { - print( - kable(mkin_long_to_wide(ds[[ds_name]]), - caption = paste("Dataset", ds_name), - booktabs = TRUE, row.names = FALSE)) - cat("\n\\clearpage\n") -} -``` - -# Parent only evaluations - -The following code performs separate fits of the candidate degradation models -to all datasets using constant variance and the two-component error model. - -```{r parent-sep, dependson = "ds"} -parent_deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB") -parent_sep_const <- mmkin( - parent_deg_mods, ds, - error_model = "const", - cluster = cl, quiet = TRUE) -parent_sep_tc <- update(parent_sep_const, error_model = "tc") -``` - -To select the parent model, the corresponding hierarchical fits are performed below. - -```{r parent-mhmkin, dependson = "parent-sep"} -parent_mhmkin <- mhmkin(list(parent_sep_const, parent_sep_tc), cluster = cl) -status(parent_mhmkin) |> kable() -``` - -All fits terminate without errors (status OK). The check for ill-defined -parameters shows that not all random effect parameters can be robustly -quantified. - -```{r dependson = "parent_mhmkin"} -illparms(parent_mhmkin) |> kable() -``` - -Therefore, the fits are updated, excluding random effects that were -ill-defined according to the `illparms` function. - -```{r parent-mhmkin-refined} -parent_mhmkin_refined <- update(parent_mhmkin, - no_random_effect = illparms(parent_mhmkin)) -status(parent_mhmkin_refined) |> kable() -``` - -The most suitable model is selected based on the AIC. - -```{r 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) -``` - -Based on the AIC, the combination of the `r best_degmod_parent` degradation -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]]) -``` - -The corresponding fit is shown below. - -```{r parent-best-full, dependson = "parent-mhmkin"} -plot(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]]) -``` - -Detailed listings of the parent fits are shown in the Appendix. - -\clearpage - -# Appendix - -## Summaries of saem fits - -### Parent fits - -```{r listings-parent, results = "asis", echo = FALSE} -for (deg_mod in parent_deg_mods) { - for (err_mod in c("const", "tc")) { - caption <- paste("Hierarchical", deg_mod, "fit with", errmods[err_mod]) - tex_listing(parent_mhmkin[[deg_mod, err_mod]], caption) - } -} -``` - -### Refined parent fits - -```{r listings-pathway, results = "asis", echo = FALSE} -for (deg_mod in parent_deg_mods) { - for (err_mod in c("const", "tc")) { - caption <- paste("Refined hierarchical", deg_mod, "fit with", errmods[err_mod]) - tex_listing(parent_mhmkin_refined[[deg_mod, err_mod]], caption) - } -} -``` - -## Session info - -```{r, echo = FALSE} -parallel::stopCluster(cl) -sessionInfo() -``` - diff --git a/inst/rmarkdown/templates/hier/template.yaml b/inst/rmarkdown/templates/hier/template.yaml deleted file mode 100644 index d8ab6a4d..00000000 --- a/inst/rmarkdown/templates/hier/template.yaml +++ /dev/null @@ -1,3 +0,0 @@ -name: Hierarchical kinetics -description: Hierarchical kinetic modelling of degradation data -create_dir: true diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex new file mode 100644 index 00000000..a2b7ce83 --- /dev/null +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex @@ -0,0 +1 @@ +\definecolor{shadecolor}{RGB}{248,248,248} diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd new file mode 100644 index 00000000..74fae164 --- /dev/null +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd @@ -0,0 +1,162 @@ +--- +title: "Hierarchical kinetic modelling of degradation data" +author: +date: +output: mkin::hierarchical_kinetics +geometry: margin=2cm +--- + +```{r packages, cache = FALSE, message = FALSE} +library(mkin) +library(knitr) +library(saemix) +library(parallel) +library(readxl) +``` + +```{r n_cores, cache = FALSE, echo = FALSE} +n_cores <- detectCores() + +if (Sys.info()["sysname"] == "Windows") { + cl <- makePSOCKcluster(n_cores) +} else { + cl <- makeForkCluster(n_cores) +} +``` + +\clearpage + +# Introduction + +This report shows hierarchical kinetic modelling for ... +The data were obtained from ... + +```{r ds} +data_path <- system.file("testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", package = "mkin") +ds <- read_spreadsheet(data_path, valid_datasets = c(1:4, 7:13)) +covariates <- attr(ds, "covariates") +``` + +The covariate data are shown below. + +```{r results = "asis", dependson = "ds", echo = FALSE} +kable(covariates, caption = "Covariate data for all datasets") +``` + +\clearpage + +The datasets with the residue time series are shown in the tables below. Please +refer to the spreadsheet for details like data sources, treatment of values +below reporting limits and time step normalisation factors. + +```{r results = "asis", dependson = "ds", echo = FALSE} +for (ds_name in names(ds)) { + print( + kable(mkin_long_to_wide(ds[[ds_name]]), + caption = paste("Dataset", ds_name), + booktabs = TRUE, row.names = FALSE)) + cat("\n\\clearpage\n") +} +``` + +# Parent only evaluations + +The following code performs separate fits of the candidate degradation models +to all datasets using constant variance and the two-component error model. + +```{r parent-sep, dependson = "ds"} +parent_deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB") +errmods <- c(const = "constant variance", tc = "two-component error") +parent_sep_const <- mmkin( + parent_deg_mods, ds, + error_model = "const", + cluster = cl, quiet = TRUE) +parent_sep_tc <- update(parent_sep_const, error_model = "tc") +``` + +To select the parent model, the corresponding hierarchical fits are performed below. + +```{r parent-mhmkin, dependson = "parent-sep"} +parent_mhmkin <- mhmkin(list(parent_sep_const, parent_sep_tc), cluster = cl) +status(parent_mhmkin) |> kable() +``` + +All fits terminate without errors (status OK). The check for ill-defined +parameters shows that not all random effect parameters can be robustly +quantified. + +```{r dependson = "parent_mhmkin"} +illparms(parent_mhmkin) |> kable() +``` + +Therefore, the fits are updated, excluding random effects that were +ill-defined according to the `illparms` function. + +```{r parent-mhmkin-refined} +parent_mhmkin_refined <- update(parent_mhmkin, + no_random_effect = illparms(parent_mhmkin)) +status(parent_mhmkin_refined) |> kable() +``` + +The most suitable model is selected based on the AIC. + +```{r, 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) +``` + +Based on the AIC, the combination of the `r best_degmod_parent` degradation +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]]) +``` + +The corresponding fit is shown below. + +```{r parent-best-full, dependson = "parent-mhmkin"} +plot(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]]) +``` + +Detailed listings of the parent fits are shown in the Appendix. + +\clearpage + +# Appendix + +## Summaries of saem fits + +### Parent fits + +```{r listings-parent, results = "asis", echo = FALSE, dependson = "parent_mhmkin"} +for (deg_mod in parent_deg_mods) { + for (err_mod in c("const", "tc")) { + caption <- paste("Hierarchical", deg_mod, "fit with", errmods[err_mod]) + tex_listing(parent_mhmkin[[deg_mod, err_mod]], caption) + } +} +``` + +### Refined parent fits + +```{r listings-parent-refined, results = "asis", echo = FALSE, dependson = "parent_mhmkin_refined"} +for (deg_mod in parent_deg_mods) { + for (err_mod in c("const", "tc")) { + caption <- paste("Refined hierarchical", deg_mod, "fit with", errmods[err_mod]) + tex_listing(parent_mhmkin_refined[[deg_mod, err_mod]], caption) + } +} +``` + +## Session info + +```{r, echo = FALSE, cache = FALSE} +parallel::stopCluster(cl) +sessionInfo() +``` + diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml b/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml new file mode 100644 index 00000000..d8ab6a4d --- /dev/null +++ b/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml @@ -0,0 +1,3 @@ +name: Hierarchical kinetics +description: Hierarchical kinetic modelling of degradation data +create_dir: true diff --git a/log/check.log b/log/check.log index 42365918..aec61e33 100644 --- a/log/check.log +++ b/log/check.log @@ -1,5 +1,5 @@ * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ -* using R version 4.2.2 (2022-10-31) +* using R version 4.2.2 Patched (2022-11-10 r83330) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 * using options ‘--no-tests --as-cran’ @@ -18,7 +18,7 @@ Maintainer: ‘Johannes Ranke ’ * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking serialization versions ... OK -* checking whether package ‘mkin’ can be installed ... [11s/11s] OK +* checking whether package ‘mkin’ can be installed ... OK * checking installed package size ... OK * checking package directory ... OK * checking for future file timestamps ... OK @@ -41,7 +41,14 @@ Maintainer: ‘Johannes Ranke ’ * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK -* checking R code for possible problems ... [19s/19s] OK +* checking R code for possible problems ... NOTE +parplot.multistart.saem.mmkin: no visible global function definition + for ‘quantile’ +Undefined global functions or variables: + quantile +Consider adding + importFrom("stats", "quantile") +to your NAMESPACE file. * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -57,7 +64,7 @@ Maintainer: ‘Johannes Ranke ’ * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... [24s/24s] OK +* checking examples ... [11s/11s] OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -69,5 +76,9 @@ Maintainer: ‘Johannes Ranke ’ * checking for detritus in the temp directory ... OK * DONE -Status: OK +Status: 1 NOTE +See + ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ +for details. + diff --git a/man/hierarchical_kinetics.Rd b/man/hierarchical_kinetics.Rd new file mode 100644 index 00000000..4bb82a4c --- /dev/null +++ b/man/hierarchical_kinetics.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hierarchical_kinetics.R +\name{hierarchical_kinetics} +\alias{hierarchical_kinetics} +\title{Hierarchical kinetics template} +\usage{ +hierarchical_kinetics(..., keep_tex = FALSE) +} +\arguments{ +\item{...}{Arguments to \code{rmarkdown::pdf_document}} + +\item{keep_tex}{Keep the intermediate tex file used in the conversion to PDF} +} +\value{ +R Markdown output format to pass to +\code{\link[rmarkdown:render]{render}} +} +\description{ +R markdown format for setting up hierarchical kinetics based on a template +provided with the mkin package. +} +\examples{ + +\dontrun{ +library(rmarkdown) +draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin") +} + +} -- cgit v1.2.3 From 9a066c13e2f6324920910a238e3b40261854e21a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 22 Dec 2022 12:57:24 +0100 Subject: Improve code formatting in template skeleton --- inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) (limited to 'inst/rmarkdown') diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd index 74fae164..3c7dc8a5 100644 --- a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd @@ -32,7 +32,9 @@ This report shows hierarchical kinetic modelling for ... The data were obtained from ... ```{r ds} -data_path <- system.file("testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", package = "mkin") +data_path <- system.file( + "testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", + package = "mkin") ds <- read_spreadsheet(data_path, valid_datasets = c(1:4, 7:13)) covariates <- attr(ds, "covariates") ``` -- cgit v1.2.3 From ed2ab7bba07e3cb036782227fe27943f4b5583fa Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 3 Jan 2023 09:31:48 +0100 Subject: Improved skeleton for hierarchical fits Now with working pathway fits using SFORB-SFO2 (only two parallel metabolites instead of three) as the data for compound Ia was not sufficient for a reliable fit. --- .gitignore | 1 + DESCRIPTION | 14 +- R/hierarchical_kinetics.R | 1 + .../hierarchical_kinetics/skeleton/skeleton.Rmd | 181 +++++++++++++++++++-- 4 files changed, 175 insertions(+), 22 deletions(-) (limited to 'inst/rmarkdown') diff --git a/.gitignore b/.gitignore index 9808896d..5416f8f0 100644 --- a/.gitignore +++ b/.gitignore @@ -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) - . Interfaces to several nonlinear - mixed-effects model packages are available, some of which are described by - Ranke et al. (2021) . Please note that no - warranty is implied for correctness of results or fitness for a particular - purpose. + . Hierarchical degradation models can + be fitted using nonlinear mixed-effects model packages as a backend as + described by Ranke et al. (2021) . 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() ``` -- cgit v1.2.3 From 18a18ef46391c4b9fa3e5e1cb06937d90a83ba05 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 3 Jan 2023 17:31:47 +0100 Subject: Improve template text, use neutral name for dll_dir --- .../templates/hierarchical_kinetics/skeleton/skeleton.Rmd | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'inst/rmarkdown') diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd index 9a975607..e26213f5 100644 --- a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd @@ -144,10 +144,10 @@ intervals are listed below. parms(parent_best, ci = TRUE) |> kable(digits = 3) ``` -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. +To investigate a potential 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, @@ -188,14 +188,14 @@ parallel formation of two metabolites is set up. ```{r m_sforb_sfo2} -if (!dir.exists("lambda_dlls")) dir.create("lambda_dlls") +if (!dir.exists("dlls")) dir.create("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", + dll_dir = "dlls", overwrite = TRUE, quiet = TRUE ) ``` -- cgit v1.2.3 From fc692554ecdaff66548cc3e6d666e44b7aaaa9af Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 3 Jan 2023 17:52:39 +0100 Subject: Neutral names for code chunks in the template --- .../templates/hierarchical_kinetics/skeleton/skeleton.Rmd | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) (limited to 'inst/rmarkdown') diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd index e26213f5..38a6bd20 100644 --- a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd +++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd @@ -186,8 +186,7 @@ parms(parent_best_pH_2, ci = TRUE) |> kable(digits = 3) 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} +```{r path-1-degmod} if (!dir.exists("dlls")) dir.create("dlls") m_sforb_sfo2 = mkinmod( @@ -203,7 +202,7 @@ m_sforb_sfo2 = mkinmod( Separate evaluations of all datasets are performed with constant variance and using two-component error. -```{r path-sep, dependson = c("m_sforb_all", "ds")} +```{r path-1-sep, dependson = c("path-1-degmod", "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") @@ -211,7 +210,7 @@ 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} +```{r dependson = "path-1-sep", fig.height = 9} plot(mixed(sforb_sep_const)) ``` @@ -219,7 +218,7 @@ 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"} +```{r path-1, dependson = "path-1-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), -- cgit v1.2.3