diff options
Diffstat (limited to 'vignettes/prebuilt/2022_dmta_pathway.rmd')
| -rw-r--r-- | vignettes/prebuilt/2022_dmta_pathway.rmd | 426 |
1 files changed, 426 insertions, 0 deletions
diff --git a/vignettes/prebuilt/2022_dmta_pathway.rmd b/vignettes/prebuilt/2022_dmta_pathway.rmd new file mode 100644 index 00000000..ff2b527c --- /dev/null +++ b/vignettes/prebuilt/2022_dmta_pathway.rmd @@ -0,0 +1,426 @@ +--- +title: "Testing hierarchical pathway kinetics with residue data on dimethenamid and dimethenamid-P" +author: Johannes Ranke +date: Last change on 8 January 2023, last compiled on `r format(Sys.time(), "%e %B %Y")` +geometry: margin=2cm +bibliography: references.bib +toc: true +output: + pdf_document: + extra_dependencies: ["float", "listing"] +--- + +```{r setup, echo = FALSE, cache = FALSE} +options(width = 80) # For summary listings +knitr::opts_chunk$set( + comment = "", tidy = FALSE, cache = TRUE, fig.pos = "H", fig.align = "center" +) +``` + +\clearpage + +# Introduction + +The purpose of this document is to test demonstrate how nonlinear hierarchical +models (NLHM) based on the parent degradation models SFO, FOMC, DFOP and HS, +with parallel formation of two or more metabolites can be fitted with the mkin package. + +It was assembled in the course of work package 1.2 of Project Number 173340 +(Application of nonlinear hierarchical models to the kinetic evaluation of +chemical degradation data) of the German Environment Agency carried out in 2022 +and 2023. + +The mkin package is used in version `r packageVersion("mkin")`, which is currently +under development. It contains the test data, and the functions used in the +evaluations. The `saemix` package is used as a backend for fitting the NLHM, +but is also loaded to make the convergence plot function available. + +This document is processed with the `knitr` package, which also provides the +`kable` function that is used to improve the display of tabular data in R +markdown documents. For parallel processing, the `parallel` package is used. + +```{r, packages, cache = FALSE, message = FALSE} +library(mkin) +library(knitr) +library(saemix) +library(parallel) +n_cores <- detectCores() +if (Sys.info()["sysname"] == "Windows") { + cl <- makePSOCKcluster(n_cores) +} else { + cl <- makeForkCluster(n_cores) +} +``` + +\clearpage + +# Data + +The test data are available in the mkin package as an object of class `mkindsg` +(mkin dataset group) under the identifier `dimethenamid_2018`. The following +preprocessing steps are done in this document. + +- The data available for the enantiomer dimethenamid-P (DMTAP) are renamed + to have the same substance name as the data for the racemic mixture + dimethenamid (DMTA). The reason for this is that no difference between their + degradation behaviour was identified in the EU risk assessment. +- Unnecessary columns are discarded +- The observation times of each dataset are multiplied with the + corresponding normalisation factor also available in the dataset, in order to + make it possible to describe all datasets with a single set of parameters + that are independent of temperature +- Finally, datasets observed in the same soil (`Elliot 1` and `Elliot 2`) are + combined, resulting in dimethenamid (DMTA) data from six soils. + +The following commented R code performs this preprocessing. + +```{r, data} +# Apply a function to each of the seven datasets in the mkindsg object to create a list +dmta_ds <- lapply(1:7, function(i) { + ds_i <- dimethenamid_2018$ds[[i]]$data # Get a dataset + ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA" # Rename DMTAP to DMTA + ds_i <- subset(ds_i, select = c("name", "time", "value")) # Select data + ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i] # Normalise time + ds_i # Return the dataset +}) + +# Use dataset titles as names for the list elements +names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title) + +# Combine data for Elliot soil to obtain a named list with six elements +dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) # +dmta_ds[["Elliot 1"]] <- NULL +dmta_ds[["Elliot 2"]] <- NULL +``` + +\clearpage + +The following tables show the `r length(dmta_ds)` datasets. + +```{r show-data, dependson = "data", results = "asis"} +for (ds_name in names(dmta_ds)) { + print( + kable(mkin_long_to_wide(dmta_ds[[ds_name]]), + caption = paste("Dataset", ds_name), + booktabs = TRUE, row.names = FALSE)) + cat("\n\\clearpage\n") +} +``` + +# Separate evaluations + +As a first step to obtain suitable starting parameters for the NLHM fits, we do +separate fits of several variants of the pathway model used previously +[@ranke2021], varying the kinetic model for the parent compound. Because the +SFORB model often provides faster convergence than the DFOP model, and can +sometimes be fitted where the DFOP model results in errors, it is included in +the set of parent models tested here. + +```{r, sep-1-const, dependson = "data"} +if (!dir.exists("dmta_dlls")) dir.create("dmta_dlls") +m_sfo_path_1 <- mkinmod( + DMTA = mkinsub("SFO", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + name = "m_sfo_path", dll_dir = "dmta_dlls", + unload = TRUE, overwrite = TRUE, + quiet = TRUE +) +m_fomc_path_1 <- mkinmod( + DMTA = mkinsub("FOMC", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + name = "m_fomc_path", dll_dir = "dmta_dlls", + unload = TRUE, overwrite = TRUE, + quiet = TRUE +) +m_dfop_path_1 <- mkinmod( + DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + name = "m_dfop_path", dll_dir = "dmta_dlls", + unload = TRUE, overwrite = TRUE, + quiet = TRUE +) +m_sforb_path_1 <- mkinmod( + DMTA = mkinsub("SFORB", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + name = "m_sforb_path", dll_dir = "dmta_dlls", + unload = TRUE, overwrite = TRUE, + quiet = TRUE +) +m_hs_path_1 <- mkinmod( + DMTA = mkinsub("HS", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + name = "m_hs_path", dll_dir = "dmta_dlls", + unload = TRUE, overwrite = TRUE, + quiet = TRUE +) +deg_mods_1 <- list( + sfo_path_1 = m_sfo_path_1, + fomc_path_1 = m_fomc_path_1, + dfop_path_1 = m_dfop_path_1, + sforb_path_1 = m_sforb_path_1, + hs_path_1 = m_hs_path_1) + +sep_1_const <- mmkin( + deg_mods_1, + dmta_ds, + error_model = "const", + quiet = TRUE) + +status(sep_1_const) |> kable() +``` + +All separate pathway fits with SFO or FOMC for the parent and constant variance +converged (status OK). Most fits with DFOP or SFORB for the parent converged +as well. The fits with HS for the parent did not converge with default settings. + +```{r, sep-1-tc, dependson = "sep-1-const"} +sep_1_tc <- update(sep_1_const, error_model = "tc") +status(sep_1_tc) |> kable() +``` + +With the two-component error model, the set of fits with convergence problems +is slightly different, with convergence problems appearing for different data +sets when applying the DFOP and SFORB model and some additional convergence +problems when using the FOMC model for the parent. + +\clearpage + +# Hierarchichal model fits + +The following code fits two sets of the corresponding hierarchical models to +the data, one assuming constant variance, and one assuming two-component error. + +```{r saem-1, dependson = c("sep-1-const", "sep-1-tc")} +saem_1 <- mhmkin(list(sep_1_const, sep_1_tc)) +``` +The run time for these fits was around two hours on five year old hardware. After +a recent hardware upgrade these fits complete in less than twenty minutes. + +```{r, saem-1-status, dependson = "saem-1"} +status(saem_1) |> kable() +``` + +According to the `status` function, all fits terminated successfully. + +```{r saem-1-anova, dependson = "saem-1"} +anova(saem_1) |> kable(digits = 1) +``` + +When the goodness-of-fit of the models is compared, a warning is obtained, +indicating that the likelihood of the pathway fit with SFORB for the parent +compound and constant variance could not be calculated with importance sampling +(method 'is'). As this is the default method on which all AIC and BIC +comparisons are based, this variant is not included in the model comparison +table. Comparing the goodness-of-fit of the remaining models, HS model model +with two-component error provides the best fit. However, for batch experiments +performed with constant conditions such as the experiments evaluated here, +there is no reason to assume a discontinuity, so the SFORB model is +preferable from a mechanistic viewpoint. In addition, the information criteria +AIC and BIC are very similar for HS and SFORB. Therefore, the SFORB model is +selected here for further refinements. + +\clearpage + +## Parameter identifiability based on the Fisher Information Matrix + +Using the `illparms` function, ill-defined statistical model parameters such as +standard deviations of the degradation parameters in the population and error +model parameters can be found. + +```{r saem-1-illparms, dependson = "saem-1"} +illparms(saem_1) |> kable() +``` + +When using constant variance, no ill-defined variance parameters are identified +with the `illparms` function in any of the degradation models. When using +the two-component error model, there is one ill-defined variance parameter +in all variants except for the variant using DFOP for the parent compound. + +For the selected combination of the SFORB pathway model with two-component +error, the random effect for the rate constant from reversibly bound DMTA to +the free DMTA (`k_DMTA_bound_free`) is not well-defined. Therefore, the fit is +updated without assuming a random effect for this parameter. + +```{r saem-sforb-path-1-tc-reduced, dependson = "saem-1"} +saem_sforb_path_1_tc_reduced <- update(saem_1[["sforb_path_1", "tc"]], + no_random_effect = "log_k_DMTA_bound_free") +illparms(saem_sforb_path_1_tc_reduced) +``` + +As expected, no ill-defined parameters remain. The model comparison below shows +that the reduced model is preferable. + +```{r saem-sforb-path-1-tc-reduced-anova, dependson = "saem-sforb-path-1-tc-reduced"} +anova(saem_1[["sforb_path_1", "tc"]], saem_sforb_path_1_tc_reduced) |> kable(digits = 1) +``` + +The convergence plot of the refined fit is shown below. + +```{r saem-sforb-path-1-tc-reduced-convergence, dependson = "saem-sforb-path-1-tc-reduced", fig.height = 12} +plot(saem_sforb_path_1_tc_reduced$so, plot.type = "convergence") +``` + +For some parameters, for example for `f_DMTA_ilr_1` and `f_DMTA_ilr_2`, i.e. +for two of the parameters determining the formation fractions of the parallel +formation of the three metabolites, some movement of the parameters is still +visible in the second phase of the algorithm. However, the amplitude of this +movement is in the range of the amplitude towards the end of the first phase. +Therefore, it is likely that an increase in iterations would not improve the +parameter estimates very much, and it is proposed that the fit is acceptable. +No numeric convergence criterion is implemented in saemix. + +\clearpage + +## Alternative check of parameter identifiability + +As an alternative check of parameter identifiability [@duchesne_2021], +multistart runs were performed on the basis of the refined fit shown above. + +```{r saem-sforb-multistart, dependson = "saem-sforb-path-1-tc-reduced"} +saem_sforb_path_1_tc_reduced_multi <- multistart(saem_sforb_path_1_tc_reduced, + n = 32, cores = 10) +``` + +```{r dependson = "saem-sforb-multistart"} +print(saem_sforb_path_1_tc_reduced_multi) +``` + +Out of the 32 fits that were initiated, only 17 terminated without an error. +The reason for this is that the wide variation of starting parameters in combination +with the parameter variation that is used in the SAEM algorithm leads to +parameter combinations for the degradation model that the numerical integration +routine cannot cope with. Because of this variation of initial parameters, +some of the model fits take up to two times more time than the original fit. + +```{r dependson = "saem-sforb-multistart", fig.cap = "Parameter boxplots for the multistart runs that succeeded", fig.height = 6, fig.width = 10} +par(mar = c(12.1, 4.1, 2.1, 2.1)) +parplot(saem_sforb_path_1_tc_reduced_multi, ylim = c(0.5, 2), las = 2) +``` + +However, visual analysis of the boxplot of the parameters obtained in the +successful fits confirms that the results are sufficiently independent of the +starting parameters, and there are no remaining ill-defined parameters. + +\clearpage + + +# Plots of selected fits + +The SFORB pathway fits with full and reduced parameter distribution model are +shown below. + +```{r fig.cap = "SFORB pathway fit with two-component error", dependson = "saem-1", fig.height = 8} +plot(saem_1[["sforb_path_1", "tc"]]) +``` + +\clearpage + +```{r fig.cap = "SFORB pathway fit with two-component error, reduced parameter model", dependson = "saem-sforb-path-1-tc-reduced", fig.height = 8} +plot(saem_sforb_path_1_tc_reduced) +``` + +Plots of the remaining fits and listings for all successful fits are shown in +the Appendix. + + +# Conclusions + +Pathway fits with SFO, FOMC, DFOP, SFORB and HS models for the parent compound +could be successfully performed. + +\clearpage + +# Acknowledgements + +The helpful comments by Janina Wöltjen of the German Environment Agency +on earlier versions of this document are gratefully acknowledged. + +# References + +\vspace{1em} + +::: {#refs} +::: + +\clearpage + +# Appendix + +## Plots of hierarchical fits not selected for refinement + +```{r fig.cap = "SFO pathway fit with two-component error", dependson = "saem-1", fig.height = 8} +plot(saem_1[["sfo_path_1", "tc"]]) +``` + +\clearpage + +```{r fig.cap = "FOMC pathway fit with two-component error", dependson = "saem-1", fig.height = 8} +plot(saem_1[["fomc_path_1", "tc"]]) +``` + +\clearpage + + +```{r fig.cap = "HS pathway fit with two-component error", dependson = "saem-1", fig.height = 8} +plot(saem_1[["sforb_path_1", "tc"]]) +``` + +\clearpage + +## Hierarchical model fit listings + +### Fits with random effects for all degradation parameters + +```{r listings-1, results = "asis", echo = FALSE} +errmods <- c(const = "constant variance", tc = "two-component error") +degmods <- c( + sfo_path_1 = "SFO path 1", + fomc_path_1 = "FOMC path 1", + dfop_path_1 = "DFOP path 1", + sforb_path_1 = "SFORB path 1", + hs_path_1 = "HS path 1") +for (deg_mod in rownames(saem_1)) { + for (err_mod in c("const", "tc")) { + fit <- saem_1[[deg_mod, err_mod]] + if (!inherits(fit$so, "try-error")) { + caption <- paste("Hierarchical", degmods[deg_mod], "fit with", errmods[err_mod]) + tex_listing(fit, caption) + } + } +} +``` + +### Improved fit of the SFORB pathway model with two-component error + +```{r listings-2, results = "asis", echo = FALSE, dependson = "listings-1"} +caption <- paste("Hierarchical SFORB pathway fit with two-component error") +tex_listing(saem_sforb_path_1_tc_reduced, caption) +``` + +## Session info + +```{r, echo = FALSE} +parallel::stopCluster(cl) +sessionInfo() +``` + +## Hardware info + +```{r, echo = FALSE} +if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) { + cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]])) +} +if(!inherits(try(meminfo <- readLines("/proc/meminfo")), "try-error")) { + cat(gsub("model name\t: ", "System memory: ", meminfo[grep("MemTotal", meminfo)[1]])) +} +``` |