From 24eb77216700cf8b2f2bde3abad84c1f83f9e32a Mon Sep 17 00:00:00 2001
From: Johannes Ranke The purpose of this document is to demonstrate how nonlinear
+hierarchical models (NLHM) based on the parent degradation models SFO,
+FOMC, DFOP and HS can be fitted with the mkin package. The mkin package is used in version 1.2.2. It contains the test data
+and the functions used in the evaluations. The This document is processed with the The test data are available in the mkin package as an object of class
+ The following commented R code performs this preprocessing. The following tables show the 6 datasets. In order to obtain suitable starting parameters for the NLHM fits,
+separate fits of the four models to the data for each soil are generated
+using the In the table above, OK indicates convergence, and C indicates failure
+to converge. All separate fits with constant variance converged, with
+the sole exception of the HS fit to the BBA 2.2 data. To prepare for
+fitting NLHM using the two-component error model, the separate fits are
+updated assuming two-component error. Using the two-component error model, the one fit that did not
+converge with constant variance did converge, but other non-SFO fits
+failed to converge. The following code fits eight versions of hierarchical models to the
+data, using SFO, FOMC, DFOP and HS for the parent compound, and using
+either constant variance or two-component error for the error model. The
+default parameter distribution model in mkin allows for variation of all
+degradation parameters across the assumed population of soils. In other
+words, each degradation parameter is associated with a random effect as
+a first step. The Convergence plots and summaries for these fits are shown in the
+appendix. The output of the The AIC and BIC values show that the biphasic models DFOP and HS give
+the best fits. The DFOP model is preferred here, as it has a better mechanistic
+basis for batch experiments with constant incubation conditions. Also,
+it shows the lowest AIC and BIC values in the first set of fits when
+combined with the two-component error model. Therefore, the DFOP model
+was selected for further refinements of the fits with the aim to make
+the model fully identifiable. Using the According to the The thus identified overparameterisation is addressed by removing the
+random effect for For the resulting fit, it is checked whether there are still
+ill-defined parameters, which is not the case. Below, the refined model is compared with the
+previous best model. The model without random effect for The AIC and BIC criteria are lower after removal of the ill-defined
+random effect for Therefore, AIC, BIC and likelihood ratio test suggest the use of the
+reduced model. The convergence of the fit is checked visually. All parameters appear to have converged to a satisfactory degree. The
+final fit is plotted using the plot method from the mkin package. Finally, a summary report of the fit is produced. The parameter check used in the The graph below shows boxplots of the parameters obtained in 50 runs
+of the saem algorithm with different parameter combinations, sampled
+from the range of the parameters obtained for the individual datasets
+fitted separately using nonlinear regression. The graph clearly confirms the lack of identifiability of the
+variance of The parameter boxplots of the multistart runs with the reduced model
+shown below indicate that all runs give similar results, regardless of
+the starting parameters. When only the parameters of the top 25% of the fits are shown (based
+on a feature introduced in mkin 1.2.2 currently under development), the
+scatter is even less as shown below. Fitting the four parent degradation models SFO, FOMC, DFOP and HS as
+part of hierarchical model fits with two different error models and
+normal distributions of the transformed degradation parameters works
+without technical problems. The biphasic models DFOP and HS gave the
+best fit to the data, but the default parameter distribution model was
+not fully identifiable. Removing the random effect for the second
+kinetic rate constant of the DFOP model resulted in a reduced model that
+was fully identifiable and showed the lowest values for the model
+selection criteria AIC and BIC. The reliability of the identification of
+all model parameters was confirmed using multiple starting values. Ranke J (2022).
+ Ranke J (2023).
mkin: Kinetic Evaluation of Chemical Degradation Data.
R package version 1.2.2, https://pkgdown.jrwb.de/mkin/.
R markdown format for setting up hierarchical kinetics based on a template
+provided with the mkin package. Arguments to Keep the intermediate tex file used in the conversion to PDF R Markdown output format to pass to
+Work package 1.1: Testing hierarchical parent
+degradation kinetics with residue data on dimethenamid and
+dimethenamid-P
+ Johannes
+Ranke
+
+ Last change on 5 January
+2022, last compiled on 5 Januar 2023
+
+ Source: vignettes/2022_wp_1.1_dmta_parent.rmd
+ 2022_wp_1.1_dmta_parent.rmd
Introduction
+
+saemix
+package is used as a backend for fitting the NLHM, but is also loaded to
+make the convergence plot function available.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.
+
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)
+}
Preprocessing of test data
+
+mkindsg
(mkin dataset group) under the identifier
+dimethenamid_2018
. The following preprocessing steps are
+still necessary:
+
+Elliot 1
+and Elliot 2
) are combined, resulting in dimethenamid
+(DMTA) data from six soils.
+
# 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, name == "DMTA", 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
+
for (ds_name in names(dmta_ds)) {
+ print(kable(mkin_long_to_wide(dmta_ds[[ds_name]]),
+ caption = paste("Dataset", ds_name),
+ label = paste0("tab:", ds_name), booktabs = TRUE))
+ cat("\n\\clearpage\n")
+}
+
+
+
+
+time
+DMTA
+
+
+0
+95.8
+
+
+0
+98.7
+
+
+14
+60.5
+
+
+30
+39.1
+
+
+59
+15.2
+
+
+120
+4.8
+
+
+
+120
+4.6
+
+
+
+
+
+time
+DMTA
+
+
+0.000000
+100.5
+
+
+0.000000
+99.6
+
+
+1.941295
+91.9
+
+
+1.941295
+91.3
+
+
+6.794534
+81.8
+
+
+6.794534
+82.1
+
+
+13.589067
+69.1
+
+
+13.589067
+68.0
+
+
+27.178135
+51.4
+
+
+27.178135
+51.4
+
+
+56.297565
+27.6
+
+
+56.297565
+26.8
+
+
+86.387643
+15.7
+
+
+86.387643
+15.3
+
+
+115.507073
+7.9
+
+
+
+115.507073
+8.1
+
+
+
+
+
+time
+DMTA
+
+
+0.0000000
+96.5
+
+
+0.0000000
+96.8
+
+
+0.0000000
+97.0
+
+
+0.6233856
+82.9
+
+
+0.6233856
+86.7
+
+
+0.6233856
+87.4
+
+
+1.8701567
+72.8
+
+
+1.8701567
+69.9
+
+
+1.8701567
+71.9
+
+
+4.3636989
+51.4
+
+
+4.3636989
+52.9
+
+
+4.3636989
+48.6
+
+
+8.7273979
+28.5
+
+
+8.7273979
+27.3
+
+
+8.7273979
+27.5
+
+
+13.0910968
+14.8
+
+
+13.0910968
+13.4
+
+
+13.0910968
+14.4
+
+
+17.4547957
+7.7
+
+
+17.4547957
+7.3
+
+
+17.4547957
+8.1
+
+
+26.1821936
+2.0
+
+
+26.1821936
+1.5
+
+
+26.1821936
+1.9
+
+
+34.9095915
+1.3
+
+
+34.9095915
+1.0
+
+
+34.9095915
+1.1
+
+
+43.6369893
+0.9
+
+
+43.6369893
+0.7
+
+
+43.6369893
+0.7
+
+
+52.3643872
+0.6
+
+
+52.3643872
+0.4
+
+
+52.3643872
+0.5
+
+
+74.8062674
+0.4
+
+
+74.8062674
+0.3
+
+
+
+74.8062674
+0.3
+
+
+
+
+
+time
+DMTA
+
+
+0.0000000
+98.09
+
+
+0.0000000
+98.77
+
+
+0.7678922
+93.52
+
+
+0.7678922
+92.03
+
+
+2.3036765
+88.39
+
+
+2.3036765
+87.18
+
+
+5.3752452
+69.38
+
+
+5.3752452
+71.06
+
+
+10.7504904
+45.21
+
+
+10.7504904
+46.81
+
+
+16.1257355
+30.54
+
+
+16.1257355
+30.07
+
+
+21.5009807
+21.60
+
+
+21.5009807
+20.41
+
+
+32.2514711
+9.10
+
+
+32.2514711
+9.70
+
+
+43.0019614
+6.58
+
+
+43.0019614
+6.31
+
+
+53.7524518
+3.47
+
+
+53.7524518
+3.52
+
+
+64.5029421
+3.40
+
+
+64.5029421
+3.67
+
+
+91.3791680
+1.62
+
+
+
+91.3791680
+1.62
+
+
+
+
+
+time
+DMTA
+
+
+0.0000000
+99.33
+
+
+0.0000000
+97.44
+
+
+0.6733938
+93.73
+
+
+0.6733938
+93.77
+
+
+2.0201814
+87.84
+
+
+2.0201814
+89.82
+
+
+4.7137565
+71.61
+
+
+4.7137565
+71.42
+
+
+9.4275131
+45.60
+
+
+9.4275131
+45.42
+
+
+14.1412696
+31.12
+
+
+14.1412696
+31.68
+
+
+18.8550262
+23.20
+
+
+18.8550262
+24.13
+
+
+28.2825393
+9.43
+
+
+28.2825393
+9.82
+
+
+37.7100523
+7.08
+
+
+37.7100523
+8.64
+
+
+47.1375654
+4.41
+
+
+47.1375654
+4.78
+
+
+56.5650785
+4.92
+
+
+56.5650785
+5.08
+
+
+80.1338612
+2.13
+
+
+
+80.1338612
+2.23
+
+
+
+
+
+time
+DMTA
+
+
+0.000000
+97.5
+
+
+0.000000
+100.7
+
+
+1.228478
+86.4
+
+
+1.228478
+88.5
+
+
+3.685435
+69.8
+
+
+3.685435
+77.1
+
+
+8.599349
+59.0
+
+
+8.599349
+54.2
+
+
+17.198697
+31.3
+
+
+17.198697
+33.5
+
+
+25.798046
+19.6
+
+
+25.798046
+20.9
+
+
+34.397395
+13.3
+
+
+34.397395
+15.8
+
+
+51.596092
+6.7
+
+
+51.596092
+8.7
+
+
+68.794789
+8.8
+
+
+68.794789
+8.7
+
+
+103.192184
+6.0
+
+
+103.192184
+4.4
+
+
+146.188928
+3.3
+
+
+146.188928
+2.8
+
+
+223.583066
+1.4
+
+
+223.583066
+1.8
+
+
+0.000000
+93.4
+
+
+0.000000
+103.2
+
+
+1.228478
+89.2
+
+
+1.228478
+86.6
+
+
+3.685435
+78.2
+
+
+3.685435
+78.1
+
+
+8.599349
+55.6
+
+
+8.599349
+53.0
+
+
+17.198697
+33.7
+
+
+17.198697
+33.2
+
+
+25.798046
+20.9
+
+
+25.798046
+19.9
+
+
+34.397395
+18.2
+
+
+34.397395
+12.7
+
+
+51.596092
+7.8
+
+
+51.596092
+9.0
+
+
+68.794789
+11.4
+
+
+68.794789
+9.0
+
+
+103.192184
+3.9
+
+
+103.192184
+4.4
+
+
+146.188928
+2.6
+
+
+146.188928
+3.4
+
+
+223.583066
+2.0
+
+
+
+223.583066
+1.7
+Separate evaluations
+
+mmkin
function from the mkin
+package. In a first step, constant variance is assumed. Convergence is
+checked with the status
function.
+
deg_mods <- c("SFO", "FOMC", "DFOP", "HS")
+f_sep_const <- mmkin(
+ deg_mods,
+ dmta_ds,
+ error_model = "const",
+ quiet = TRUE)
+
+status(f_sep_const) |> kable()
+
+
+
+
+
+ Calke
+Borstel
+Flaach
+BBA 2.2
+BBA 2.3
+Elliot
+
+
+SFO
+OK
+OK
+OK
+OK
+OK
+OK
+
+
+FOMC
+OK
+OK
+OK
+OK
+OK
+OK
+
+
+DFOP
+OK
+OK
+OK
+OK
+OK
+OK
+
+
+
+HS
+OK
+OK
+OK
+C
+OK
+OK
+
+
+
+
+
+
+ Calke
+Borstel
+Flaach
+BBA 2.2
+BBA 2.3
+Elliot
+
+
+SFO
+OK
+OK
+OK
+OK
+OK
+OK
+
+
+FOMC
+OK
+OK
+OK
+OK
+C
+OK
+
+
+DFOP
+OK
+OK
+C
+OK
+C
+OK
+
+
+
+HS
+OK
+C
+OK
+OK
+OK
+OK
+Hierarchichal model fits
+
+mhmkin
function makes it possible to fit
+all eight versions in parallel (given a sufficient number of computing
+cores being available) to save execution time.status
function shows that all fits
+terminated successfully.
+
+
+
+
+
+ const
+tc
+
+
+SFO
+OK
+OK
+
+
+FOMC
+OK
+OK
+
+
+DFOP
+OK
+OK
+
+
+
+HS
+OK
+OK
+
+
+
+
+
+
+ npar
+AIC
+BIC
+Lik
+
+
+SFO const
+5
+796.3
+795.3
+-393.2
+
+
+SFO tc
+6
+798.3
+797.1
+-393.2
+
+
+FOMC const
+7
+734.2
+732.7
+-360.1
+
+
+FOMC tc
+8
+720.4
+718.8
+-352.2
+
+
+DFOP const
+9
+711.8
+710.0
+-346.9
+
+
+HS const
+9
+714.0
+712.1
+-348.0
+
+
+DFOP tc
+10
+665.5
+663.4
+-322.8
+
+
+
+HS tc
+10
+667.1
+665.0
+-323.6
+Parameter identifiability based on the Fisher Information
+Matrix
+
+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.
+
+
+
+
+
+ const
+tc
+
+
+SFO
+
+ b.1
+
+
+FOMC
+
+ sd(DMTA_0)
+
+
+DFOP
+sd(k2)
+sd(k2)
+
+
+
+HS
+
+ sd(tb)
+illparms
function, the fitted standard
+deviation of the second kinetic rate constant k2
is
+ill-defined in both DFOP fits. This suggests that different values would
+be obtained for this standard deviation when using different starting
+values.k2
from the parameter model.
+
f_saem_dfop_tc_no_ranef_k2 <- update(f_saem[["DFOP", "tc"]],
+ no_random_effect = "k2")
+
illparms(f_saem_dfop_tc_no_ranef_k2)
k2
+is a reduced version of the previous model. Therefore, the models are
+nested and can be compared using the likelihood ratio test. This is
+achieved with the argument test = TRUE
to the
+anova
function.
+
anova(f_saem[["DFOP", "tc"]], f_saem_dfop_tc_no_ranef_k2, test = TRUE) |>
+ kable(format.args = list(digits = 4))
+
+
+
+
+
+ npar
+AIC
+BIC
+Lik
+Chisq
+Df
+Pr(>Chisq)
+
+
+f_saem_dfop_tc_no_ranef_k2
+9
+663.8
+661.9
+-322.9
+NA
+NA
+NA
+
+
+
+f_saem[[“DFOP”, “tc”]]
+10
+665.5
+663.4
+-322.8
+0.2809
+1
+0.5961
+k2
. The p value of the likelihood ratio
+test is much greater than 0.05, indicating that the model with the
+higher likelihood (here the model with random effects for all
+degradation parameters f_saem[["DFOP", "tc"]]
) does not fit
+significantly better than the model with the lower likelihood (the
+reduced model f_saem_dfop_tc_no_ranef_k2
).
+
plot(f_saem_dfop_tc_no_ranef_k2)
+
summary(f_saem_dfop_tc_no_ranef_k2)
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:13 2023
+Date of summary: Thu Jan 5 08:19:13 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 4.075 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.759266 0.087034 0.009933 0.930827
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.76 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 663.8 661.9 -322.9
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.228939 96.285869 100.17201
+k1 0.064063 0.033477 0.09465
+k2 0.008297 0.005824 0.01077
+g 0.953821 0.914328 0.99331
+a.1 1.068479 0.869538 1.26742
+b.1 0.029424 0.022406 0.03644
+SD.DMTA_0 2.030437 0.404824 3.65605
+SD.k1 0.594692 0.256660 0.93272
+SD.g 1.006754 0.361327 1.65218
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0218
+k2 0.0556 0.0355
+g -0.0516 -0.0284 -0.2800
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.0304 0.4048 3.6560
+SD.k1 0.5947 0.2567 0.9327
+SD.g 1.0068 0.3613 1.6522
+
+Variance model:
+ est. lower upper
+a.1 1.06848 0.86954 1.26742
+b.1 0.02942 0.02241 0.03644
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.45 41.4 12.46 10.82 83.54
Alternative check of parameter identifiability
+
+illparms
function is
+based on a quadratic approximation of the likelihood surface near its
+optimum, which is calculated using the Fisher Information Matrix (FIM).
+An alternative way to check parameter identifiability based on a
+multistart approach has recently been implemented in mkin.
+
f_saem_dfop_tc_multi <- multistart(f_saem[["DFOP", "tc"]], n = 50, cores = 15)
+
par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_multi, lpos = "bottomright", ylim = c(0.3, 10), las = 2)
k2
in the full model. The overparameterisation
+of the model also indicates a lack of identifiability of the variance of
+parameter g
.
+
f_saem_dfop_tc_no_ranef_k2_multi <- multistart(f_saem_dfop_tc_no_ranef_k2,
+ n = 50, cores = 15)
+
par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2,
+ lpos = "bottomright")
+
par(mar = c(6.1, 4.1, 2.1, 2.1))
+parplot(f_saem_dfop_tc_no_ranef_k2_multi, ylim = c(0.5, 2), las = 2, llquant = 0.25,
+ lpos = "bottomright")
Conclusions
+
+Appendix
+
+Hierarchical model fit listings
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:06 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - k_DMTA * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 1.09 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k_DMTA
+97.2953 0.0566
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k_DMTA
+DMTA_0 97.3 0
+k_DMTA 0.0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 796.3 795.3 -393.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.28130 95.71113 98.8515
+k_DMTA 0.05665 0.02909 0.0842
+a.1 2.66442 2.35579 2.9731
+SD.DMTA_0 1.54776 0.15447 2.9411
+SD.k_DMTA 0.60690 0.26248 0.9513
+
+Correlation:
+ DMTA_0
+k_DMTA 0.0168
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5478 0.1545 2.9411
+SD.k_DMTA 0.6069 0.2625 0.9513
+
+Variance model:
+ est. lower upper
+a.1 2.664 2.356 2.973
+
+Estimated disappearance times:
+ DT50 DT90
+DMTA 12.24 40.65
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - k_DMTA * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.441 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k_DMTA
+96.99175 0.05603
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k_DMTA
+DMTA_0 96.99 0
+k_DMTA 0.00 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 798.3 797.1 -393.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.271822 95.703157 98.84049
+k_DMTA 0.056638 0.029110 0.08417
+a.1 2.660081 2.230398 3.08976
+b.1 0.001665 -0.006911 0.01024
+SD.DMTA_0 1.545520 0.145035 2.94601
+SD.k_DMTA 0.606422 0.262274 0.95057
+
+Correlation:
+ DMTA_0
+k_DMTA 0.0169
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5455 0.1450 2.9460
+SD.k_DMTA 0.6064 0.2623 0.9506
+
+Variance model:
+ est. lower upper
+a.1 2.660081 2.230398 3.08976
+b.1 0.001665 -0.006911 0.01024
+
+Estimated disappearance times:
+ DT50 DT90
+DMTA 12.24 40.65
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:06 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 1.156 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 alpha beta
+ 98.292 9.909 156.341
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 alpha beta
+DMTA_0 98.29 0 0
+alpha 0.00 1 0
+beta 0.00 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 734.2 732.7 -360.1
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.3435 96.9033 99.784
+alpha 7.2007 2.5889 11.812
+beta 112.8746 34.8816 190.868
+a.1 2.0459 1.8054 2.286
+SD.DMTA_0 1.4795 0.2717 2.687
+SD.alpha 0.6396 0.1509 1.128
+SD.beta 0.6874 0.1587 1.216
+
+Correlation:
+ DMTA_0 alpha
+alpha -0.1125
+beta -0.1227 0.3632
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.4795 0.2717 2.687
+SD.alpha 0.6396 0.1509 1.128
+SD.beta 0.6874 0.1587 1.216
+
+Variance model:
+ est. lower upper
+a.1 2.046 1.805 2.286
+
+Estimated disappearance times:
+ DT50 DT90 DT50back
+DMTA 11.41 42.53 12.8
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.729 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+DMTA_0 alpha beta
+98.772 4.663 92.597
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 alpha beta
+DMTA_0 98.77 0 0
+alpha 0.00 1 0
+beta 0.00 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 720.4 718.8 -352.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.99136 97.26011 100.72261
+alpha 5.86312 2.57485 9.15138
+beta 88.55571 29.20889 147.90254
+a.1 1.51063 1.24384 1.77741
+b.1 0.02824 0.02040 0.03609
+SD.DMTA_0 1.57436 -0.04867 3.19739
+SD.alpha 0.59871 0.17132 1.02611
+SD.beta 0.72994 0.22849 1.23139
+
+Correlation:
+ DMTA_0 alpha
+alpha -0.1363
+beta -0.1414 0.2542
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5744 -0.04867 3.197
+SD.alpha 0.5987 0.17132 1.026
+SD.beta 0.7299 0.22849 1.231
+
+Variance model:
+ est. lower upper
+a.1 1.51063 1.2438 1.77741
+b.1 0.02824 0.0204 0.03609
+
+Estimated disappearance times:
+ DT50 DT90 DT50back
+DMTA 11.11 42.6 12.82
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.007 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.64383 0.09211 0.02999 0.76814
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.64 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 711.8 710 -346.9
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.092481 96.573898 99.61106
+k1 0.062499 0.030336 0.09466
+k2 0.009065 -0.005133 0.02326
+g 0.948967 0.862079 1.03586
+a.1 1.821671 1.604774 2.03857
+SD.DMTA_0 1.677785 0.472066 2.88350
+SD.k1 0.634962 0.270788 0.99914
+SD.k2 1.033498 -0.205994 2.27299
+SD.g 1.710046 0.428642 2.99145
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0246
+k2 0.0491 0.0953
+g -0.0552 -0.0889 -0.4795
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.678 0.4721 2.8835
+SD.k1 0.635 0.2708 0.9991
+SD.k2 1.033 -0.2060 2.2730
+SD.g 1.710 0.4286 2.9914
+
+Variance model:
+ est. lower upper
+a.1 1.822 1.605 2.039
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.79 42.8 12.88 11.09 76.46
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:08 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 3.033 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.759266 0.087034 0.009933 0.930827
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.76 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 665.5 663.4 -322.8
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.377019 96.447952 100.30609
+k1 0.064843 0.034607 0.09508
+k2 0.008895 0.006368 0.01142
+g 0.949696 0.903815 0.99558
+a.1 1.065241 0.865754 1.26473
+b.1 0.029340 0.022336 0.03634
+SD.DMTA_0 2.007754 0.387982 3.62753
+SD.k1 0.580473 0.250286 0.91066
+SD.k2 0.006105 -4.920337 4.93255
+SD.g 1.097149 0.412779 1.78152
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0235
+k2 0.0595 0.0424
+g -0.0470 -0.0278 -0.2731
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.007754 0.3880 3.6275
+SD.k1 0.580473 0.2503 0.9107
+SD.k2 0.006105 -4.9203 4.9325
+SD.g 1.097149 0.4128 1.7815
+
+Variance model:
+ est. lower upper
+a.1 1.06524 0.86575 1.26473
+b.1 0.02934 0.02234 0.03634
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.36 41.32 12.44 10.69 77.92
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.004 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 tb
+97.82176 0.06931 0.02997 11.13945
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 tb
+DMTA_0 97.82 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+tb 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 714 712.1 -348
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.16102 96.47747 99.84456
+k1 0.07876 0.05261 0.10491
+k2 0.02227 0.01706 0.02747
+tb 13.99089 -7.40049 35.38228
+a.1 1.82305 1.60700 2.03910
+SD.DMTA_0 1.88413 0.56204 3.20622
+SD.k1 0.34292 0.10482 0.58102
+SD.k2 0.19851 0.01718 0.37985
+SD.tb 1.68168 0.58064 2.78272
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0142
+k2 0.0001 -0.0025
+tb 0.0165 -0.1256 -0.0301
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.8841 0.56204 3.2062
+SD.k1 0.3429 0.10482 0.5810
+SD.k2 0.1985 0.01718 0.3798
+SD.tb 1.6817 0.58064 2.7827
+
+Variance model:
+ est. lower upper
+a.1 1.823 1.607 2.039
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 8.801 67.91 20.44 8.801 31.13
+
+
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:08 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 3.287 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 tb
+98.45190 0.07525 0.02576 19.19375
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 tb
+DMTA_0 98.45 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+tb 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 667.1 665 -323.6
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.76570 95.81350 99.71791
+k1 0.05855 0.03080 0.08630
+k2 0.02337 0.01664 0.03010
+tb 31.09638 29.38289 32.80987
+a.1 1.08835 0.88590 1.29080
+b.1 0.02964 0.02257 0.03671
+SD.DMTA_0 2.04877 0.42607 3.67147
+SD.k1 0.59166 0.25621 0.92711
+SD.k2 0.30698 0.09561 0.51835
+SD.tb 0.01274 -0.10914 0.13462
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0160
+k2 -0.0070 -0.0024
+tb -0.0668 -0.0103 -0.2013
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.04877 0.42607 3.6715
+SD.k1 0.59166 0.25621 0.9271
+SD.k2 0.30698 0.09561 0.5183
+SD.tb 0.01274 -0.10914 0.1346
+
+Variance model:
+ est. lower upper
+a.1 1.08835 0.88590 1.29080
+b.1 0.02964 0.02257 0.03671
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.84 51.71 15.57 11.84 29.66
+
+
Session info
+
+
+R version 4.2.2 Patched (2022-11-10 r83330)
+Platform: x86_64-pc-linux-gnu (64-bit)
+Running under: Debian GNU/Linux bookworm/sid
+
+Matrix products: default
+BLAS: /usr/lib/x86_64-linux-gnu/openblas-serial/libblas.so.3
+LAPACK: /usr/lib/x86_64-linux-gnu/openblas-serial/libopenblas-r0.3.21.so
+
+locale:
+ [1] LC_CTYPE=de_DE.UTF-8 LC_NUMERIC=C
+ [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=de_DE.UTF-8
+ [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=de_DE.UTF-8
+ [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
+ [9] LC_ADDRESS=C LC_TELEPHONE=C
+[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
+
+attached base packages:
+[1] parallel stats graphics grDevices utils datasets methods
+[8] base
+
+other attached packages:
+[1] saemix_3.2 npde_3.3 knitr_1.41 mkin_1.2.2
+
+loaded via a namespace (and not attached):
+ [1] deSolve_1.34 zoo_1.8-11 tidyselect_1.2.0 xfun_0.35
+ [5] bslib_0.4.2 purrr_1.0.0 lattice_0.20-45 colorspace_2.0-3
+ [9] vctrs_0.5.1 generics_0.1.3 htmltools_0.5.4 yaml_2.3.6
+[13] utf8_1.2.2 rlang_1.0.6 pkgdown_2.0.7 jquerylib_0.1.4
+[17] pillar_1.8.1 glue_1.6.2 DBI_1.1.3 lifecycle_1.0.3
+[21] stringr_1.5.0 munsell_0.5.0 gtable_0.3.1 ragg_1.2.4
+[25] codetools_0.2-18 memoise_2.0.1 evaluate_0.19 fastmap_1.1.0
+[29] lmtest_0.9-40 fansi_1.0.3 highr_0.9 scales_1.2.1
+[33] cachem_1.0.6 desc_1.4.2 jsonlite_1.8.4 systemfonts_1.0.4
+[37] fs_1.5.2 textshaping_0.3.6 gridExtra_2.3 ggplot2_3.4.0
+[41] digest_0.6.31 stringi_1.7.8 dplyr_1.0.10 grid_4.2.2
+[45] rprojroot_2.0.3 cli_3.5.0 tools_4.2.2 magrittr_2.0.3
+[49] sass_0.4.4 tibble_3.1.8 pkgconfig_2.0.3 assertthat_0.2.1
+[53] rmarkdown_2.19 R6_2.5.1 mclust_6.0.0 nlme_3.1-161
+[57] compiler_4.2.2
All vignettes
- @Manual{,
title = {mkin: Kinetic Evaluation of Chemical Degradation Data},
author = {Johannes Ranke},
- year = {2022},
+ year = {2023},
note = {R package version 1.2.2},
url = {https://pkgdown.jrwb.de/mkin/},
}
@@ -136,7 +136,7 @@ R package version 1.2.2, https://pkgdown
diff --git a/docs/dev/index.html b/docs/dev/index.html
index 4723879e..993b8eea 100644
--- a/docs/dev/index.html
+++ b/docs/dev/index.html
@@ -19,11 +19,11 @@
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.">
hierarchical_kinetics(..., keep_tex = FALSE)
Arguments
+ rmarkdown::pdf_document
Value
+
+
+render
Create and work with nonlinear hierarchical models
Hierarchical kinetics template
Read datasets and relevant meta information from a spreadsheet file
Wrap the output of a summary function in tex listing environment
Display the output of a summary function according to the output format
mkinsub()
has an argument to
, specifying names of
variables to which a transfer is to be assumed in the model.
If the argument use_of_ff
is set to "min"
-(default) and the model for the compartment is "SFO" or "SFORB", an
+and the model for the compartment is "SFO" or "SFORB", an
additional mkinsub()
argument can be sink = FALSE
, effectively
fixing the flux to sink to zero.
In print.mkinmod, this argument is currently not used.
@@ -247,7 +247,7 @@ in the FOCUS and NAFTA guidance documents are used.
For kinetic models with more than one observed variable, a symbolic solution of the system of differential equations is included in the resulting mkinmod object in some cases, speeding up the solution.
-If a C compiler is found by pkgbuild::has_compiler()
and there
+
If a C compiler is found by pkgbuild::has_compiler()
and there
is more than one observed variable in the specification, C code is generated
for evaluating the differential equations, compiled using
inline::cfunction()
and added to the resulting mkinmod object.
R/summary_listing.R
+ summary_listing.Rd
This function is intended for use in a R markdown code chunk with the chunk
+option results = "asis"
.
summary_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+tex_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+html_listing(object, caption = NULL)
The object for which the summary is to be listed
An optional caption
An optional label, ignored in html output
Should a new page be started after the listing? Ignored in html output