From 886c9ef013124aa954d960c655b349b5340ff154 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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.
---
 .../hierarchical_kinetics/skeleton/skeleton.Rmd    | 162 +++++++++++++++++++++
 1 file changed, 162 insertions(+)
 create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd

(limited to 'inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd')

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()
+```
+
-- 
cgit v1.2.1


From 9a066c13e2f6324920910a238e3b40261854e21a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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/templates/hierarchical_kinetics/skeleton/skeleton.Rmd')

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.1


From ed2ab7bba07e3cb036782227fe27943f4b5583fa Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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.
---
 .../hierarchical_kinetics/skeleton/skeleton.Rmd    | 181 +++++++++++++++++++--
 1 file changed, 166 insertions(+), 15 deletions(-)

(limited to 'inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd')

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.1


From 18a18ef46391c4b9fa3e5e1cb06937d90a83ba05 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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/templates/hierarchical_kinetics/skeleton/skeleton.Rmd')

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.1


From fc692554ecdaff66548cc3e6d666e44b7aaaa9af Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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/templates/hierarchical_kinetics/skeleton/skeleton.Rmd')

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.1

