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.

4 files changed, 175 insertions, 22 deletions

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