From 5364f037a72863ef5ba81e14ba4417f68fd389f9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 18 Nov 2022 19:14:47 +0100
Subject: Make mixed model test data permanent to ensure reproducibility

To ensure that tests on different platforms work on the same data, the
mixed modelling test data previosly generated in
tests/testthat/setup_script.R were generated once using the script in
inst/dataset/generation/ds_mixed.R, and are now distributed with the
package.
---
 inst/dataset_generation/ds_mixed.R | 105 +++++++++++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)
 create mode 100644 inst/dataset_generation/ds_mixed.R

(limited to 'inst')

diff --git a/inst/dataset_generation/ds_mixed.R b/inst/dataset_generation/ds_mixed.R
new file mode 100644
index 00000000..f2ae6e7e
--- /dev/null
+++ b/inst/dataset_generation/ds_mixed.R
@@ -0,0 +1,105 @@
+# Synthetic data for hierarchical kinetic models
+# Refactored version of the code previously in tests/testthat/setup_script.R
+# The number of datasets was 3 for FOMC, and 10 for HS in that script, now it
+# is always 15 for consistency
+
+library(mkin)  # We use mkinmod and mkinpredict
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+n <- 15
+log_sd <- 0.3
+err_1 = list(const = 1, prop = 0.05)
+tc <- function(value) sigma_twocomp(value, err_1$const, err_1$prop)
+const <- function(value) 2
+
+set.seed(123456)
+SFO <- mkinmod(parent = mkinsub("SFO"))
+sfo_pop <- list(parent_0 = 100, k_parent = 0.03)
+sfo_parms <- as.matrix(data.frame(
+    k_parent = rlnorm(n, log(sfo_pop$k_parent), log_sd)))
+set.seed(123456)
+ds_sfo <- lapply(1:n, function(i) {
+  ds_mean <- mkinpredict(SFO, sfo_parms[i, ],
+    c(parent = sfo_pop$parent_0), sampling_times)
+  add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_sfo, "pop") <- sfo_pop
+attr(ds_sfo, "parms") <- sfo_parms
+
+set.seed(123456)
+FOMC <- mkinmod(parent = mkinsub("FOMC"))
+fomc_pop <- list(parent_0 = 100, alpha = 2, beta = 8)
+fomc_parms <- as.matrix(data.frame(
+    alpha = rlnorm(n, log(fomc_pop$alpha), 0.4),
+    beta = rlnorm(n, log(fomc_pop$beta), 0.2)))
+set.seed(123456)
+ds_fomc <- lapply(1:n, function(i) {
+  ds_mean <- mkinpredict(FOMC, fomc_parms[i, ],
+    c(parent = fomc_pop$parent_0), sampling_times)
+  add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_fomc, "pop") <- fomc_pop
+attr(ds_fomc, "parms") <- fomc_parms
+
+set.seed(123456)
+DFOP <- mkinmod(parent = mkinsub("DFOP"))
+dfop_pop <- list(parent_0 = 100, k1 = 0.06, k2 = 0.015, g = 0.4)
+dfop_parms <- as.matrix(data.frame(
+  k1 = rlnorm(n, log(dfop_pop$k1), log_sd),
+  k2 = rlnorm(n, log(dfop_pop$k2), log_sd),
+  g = plogis(rnorm(n, qlogis(dfop_pop$g), log_sd))))
+set.seed(123456)
+ds_dfop <- lapply(1:n, function(i) {
+  ds_mean <- mkinpredict(DFOP, dfop_parms[i, ],
+    c(parent = dfop_pop$parent_0), sampling_times)
+  add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_dfop, "pop") <- dfop_pop
+attr(ds_dfop, "parms") <- dfop_parms
+
+set.seed(123456)
+HS <- mkinmod(parent = mkinsub("HS"))
+hs_pop <- list(parent_0 = 100, k1 = 0.08, k2 = 0.01, tb = 15)
+hs_parms <- as.matrix(data.frame(
+  k1 = rlnorm(n, log(hs_pop$k1), log_sd),
+  k2 = rlnorm(n, log(hs_pop$k2), log_sd),
+  tb = rlnorm(n, log(hs_pop$tb), 0.1)))
+set.seed(123456)
+ds_hs <- lapply(1:n, function(i) {
+  ds_mean <- mkinpredict(HS, hs_parms[i, ],
+    c(parent = hs_pop$parent_0), sampling_times)
+  add_err(ds_mean, const, n = 1)[[1]]
+})
+attr(ds_hs, "pop") <- hs_pop
+attr(ds_hs, "parms") <- hs_parms
+
+set.seed(123456)
+DFOP_SFO <- mkinmod(
+  parent = mkinsub("DFOP", "m1"),
+  m1 = mkinsub("SFO"),
+  quiet = TRUE)
+dfop_sfo_pop <- list(parent_0 = 100,
+  k_m1 = 0.007, f_parent_to_m1 = 0.5,
+  k1 = 0.1, k2 = 0.02, g = 0.5)
+dfop_sfo_parms <- as.matrix(data.frame(
+  k1 = rlnorm(n, log(dfop_sfo_pop$k1), log_sd),
+  k2 = rlnorm(n, log(dfop_sfo_pop$k2), log_sd),
+  g = plogis(rnorm(n, qlogis(dfop_sfo_pop$g), log_sd)),
+  f_parent_to_m1 = plogis(rnorm(n,
+      qlogis(dfop_sfo_pop$f_parent_to_m1), log_sd)),
+  k_m1 = rlnorm(n, log(dfop_sfo_pop$k_m1), log_sd)))
+ds_dfop_sfo_mean <- lapply(1:n,
+  function(i) {
+    mkinpredict(DFOP_SFO, dfop_sfo_parms[i, ],
+      c(parent = dfop_sfo_pop$parent_0, m1 = 0), sampling_times)
+  }
+)
+set.seed(123456)
+ds_dfop_sfo <- lapply(ds_dfop_sfo_mean, function(ds) {
+  add_err(ds,
+    sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2),
+    n = 1, secondary = "m1")[[1]]
+})
+attr(ds_dfop_sfo, "pop") <- dfop_sfo_pop
+attr(ds_dfop_sfo, "parms") <- dfop_sfo_parms
+
+#save(ds_sfo, ds_fomc, ds_dfop, ds_hs, ds_dfop_sfo, file = "data/ds_mixed.rda", version = 2)
-- 
cgit v1.2.1


From 0023df3c31fac29b5f9337ecd732a5dfd4d51a2d Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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 +
 .../lambda-cyhalothrin_soil_efsa_2014.xlsx         | Bin 0 -> 36231 bytes
 3 files changed, 189 insertions(+)
 create mode 100644 inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd
 create mode 100644 inst/rmarkdown/templates/hier/template.yaml
 create mode 100644 inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx

(limited to 'inst')

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
diff --git a/inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx b/inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx
new file mode 100644
index 00000000..32fc049f
Binary files /dev/null and b/inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx differ
-- 
cgit v1.2.1


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.
---
 .../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 +
 5 files changed, 166 insertions(+), 189 deletions(-)
 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

(limited to 'inst')

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
-- 
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')

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')

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')

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')

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


From 24eb77216700cf8b2f2bde3abad84c1f83f9e32a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 9 Jan 2023 06:22:04 +0100
Subject: Prebuilt PDF vignettes, summary_listing

---
 inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx | Bin 0 -> 35878 bytes
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx

(limited to 'inst')

diff --git a/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx b/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx
new file mode 100644
index 00000000..3252fdf1
Binary files /dev/null and b/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx differ
-- 
cgit v1.2.1