\clearpage
library(mkin) library(knitr) library(saemix) library(parallel)
n_cores <- detectCores() # We need to start a new cluster after defining a compiled model that is # saved as a DLL to the user directory, therefore we define a function # This is used again after defining the pathway model start_cluster <- function(n_cores) { if (Sys.info()["sysname"] == "Windows") { ret <- makePSOCKcluster(n_cores) } else { ret <- makeForkCluster(n_cores) } return(ret) } cl <- start_cluster(n_cores)
\clearpage
This report shows hierarchical kinetic modelling for ... The data were obtained from ...
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.
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.
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") }
The following code performs separate fits of the candidate degradation models to all datasets using constant variance and the two-component error model.
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.
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.
illparms(parent_mhmkin) |> kable()
Therefore, the fits are updated, excluding random effects that were
ill-defined according to the illparms
function. The status of the fits
is checked.
parent_mhmkin_refined <- update(parent_mhmkin, no_random_effect = illparms(parent_mhmkin)) status(parent_mhmkin_refined) |> kable()
Also, it is checked if the AIC values of the refined fits are actually smaller than the AIC values of the original fits.
(AIC(parent_mhmkin_refined) < AIC(parent_mhmkin)) |> kable()
From the refined fits, the most suitable model is selected using the AIC.
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
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.
illparms(parent_best)
The corresponding fit is plotted below.
plot(parent_best)
The fitted parameters, together with approximate confidence intervals are listed below.
parms(parent_best, ci = TRUE) |> kable(digits = 3)
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.
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.
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.
plot(parent_best_pH_2)
illparms(parent_best_pH_2) parms(parent_best_pH_2, ci = TRUE) |> kable(digits = 3)
The endpoints corresponding to the median pH in the tested soils are shown below.
endpoints(parent_best_pH_2)
stopCluster(cl)
\clearpage
As an example of a pathway fit, a model with SFORB for the parent compound and parallel formation of two metabolites is set up.
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 = "dlls", overwrite = TRUE, quiet = TRUE )
Separate evaluations of all datasets are performed with constant variance and using two-component error.
cl_path <- NULL # Uncomment the next line for parallel processing. This does not work on Windows, # as stored dlls (that we need for caching to work) do not seem to work on # PSOCK clusters generated on Windows. # cl_path <- start_cluster(n_cores) sforb_sep_const <- mmkin(list(sforb_path = m_sforb_sfo2), ds, cluster = cl_path, quiet = TRUE) sforb_sep_tc <- update(sforb_sep_const, error_model = "tc")
The separate fits with constant variance are plotted.
plot(mixed(sforb_sep_const))
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.
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_path)
status(path_1) |> kable()
The status information shows that both fits were successfully completed.
anova(path_1) |> kable(digits = 1)
Model comparison shows that the two-component error model provides a much better fit.
illparms(path_1[["sforb_path", "tc"]])
Two ill-defined variance components are found. Therefore, the fit is repeated with the corresponding random effects removed.
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.
illparms(path_1_refined)
Below, the refined fit is plotted and the fitted parameters are shown together with their 95% confidence intervals.
plot(path_1_refined)
parms(path_1_refined, ci = TRUE) |> kable(digits = 3)
The pathway endpoints corresponding to the median pH in the tested soils are shown below.
endpoints(path_1_refined)
if (!is.null(cl_path)) stopCluster(cl_path)
\clearpage
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) } }
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) } }
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")
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.