Introduction

+ + +

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 can be fitted with the mkin package, also considering @@ -377,7 +103,7 @@ parameters. Because in some other case studies, the SFORB parameterisation of biexponential decline has shown some advantages over the DFOP parameterisation, SFORB was included in the list of tested models as well.

The mkin package is used in version 1.2.9, which is contains the +

The mkin package is used in version 1.2.10, which is contains the functions that were used for 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.

@@ -385,33 +111,35 @@ make the convergence plot function available.

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

Test data

data_file <- system.file(
-  "testdata", "mesotrione_soil_efsa_2016.xlsx", package = "mkin")
-meso_ds <- read_spreadsheet(data_file, parent_only = TRUE)

+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)
+}

Test data +

+data_file <- system.file(
+  "testdata", "mesotrione_soil_efsa_2016.xlsx", package = "mkin")
+meso_ds <- read_spreadsheet(data_file, parent_only = TRUE)

The following tables show the covariate data and the 18 datasets that were read in from the spreadsheet file.

pH <- attr(meso_ds, "covariates")
-kable(pH, caption = "Covariate data")

- +

+pH <- attr(meso_ds, "covariates")
+kable(pH, caption = "Covariate data")

- - + - - + @@ -487,20 +215,19 @@ kable(pH, caption = "Covariate data")

Covariate data
	pH
Richmond

for (ds_name in names(meso_ds)) {
-  print(
-    kable(mkin_long_to_wide(meso_ds[[ds_name]]),
-      caption = paste("Dataset", ds_name),
-      booktabs = TRUE, row.names = FALSE))
-}

- +

+for (ds_name in names(meso_ds)) {
+  print(
+    kable(mkin_long_to_wide(meso_ds[[ds_name]]),
+      caption = paste("Dataset", ds_name),
+      booktabs = TRUE, row.names = FALSE))
+}

- - + - - + @@ -592,14 +319,12 @@ kable(pH, caption = "Covariate data")

Dataset Richmond
time	meso
0.000000

- +

- - + - - + @@ -635,14 +360,12 @@ kable(pH, caption = "Covariate data")

Dataset Richmond 2
time	meso
0.000000

- +

- - + - - + @@ -678,14 +401,12 @@ kable(pH, caption = "Covariate data")

Dataset ERTC
time	meso
0.000000

- +

- - + - - + @@ -717,14 +438,12 @@ kable(pH, caption = "Covariate data")

Dataset Toulouse
time	meso
0.000000

- +

- - + - - + @@ -752,14 +471,12 @@ kable(pH, caption = "Covariate data")

Dataset Picket Piece
time	meso
0.000000

- +

- - + - - + @@ -783,14 +500,12 @@ kable(pH, caption = "Covariate data")

Dataset 721
time	meso
0.00000

- +

- - + - - + @@ -814,14 +529,12 @@ kable(pH, caption = "Covariate data")

Dataset 722
time	meso
0.00000

- +

- - + - - + @@ -845,14 +558,12 @@ kable(pH, caption = "Covariate data")

Dataset 723
time	meso
0.00000

- +

- - + - - + @@ -876,14 +587,12 @@ kable(pH, caption = "Covariate data")

Dataset 724
time	meso
0.000000

- +

- - + - - + @@ -907,14 +616,12 @@ kable(pH, caption = "Covariate data")

Dataset 725
time	meso
0.00000

- +

- - + - - + @@ -938,14 +645,12 @@ kable(pH, caption = "Covariate data")

Dataset 727
time	meso
0.00000

- +

- - + - - + @@ -969,14 +674,12 @@ kable(pH, caption = "Covariate data")

Dataset 728
time	meso
0.00000

- +

- - + - - + @@ -1000,14 +703,12 @@ kable(pH, caption = "Covariate data")

Dataset 729
time	meso
0.00000

- +

- - + - - + @@ -1031,14 +732,12 @@ kable(pH, caption = "Covariate data")

Dataset 730
time	meso
0.00000

- +

- - + - - + @@ -1062,14 +761,12 @@ kable(pH, caption = "Covariate data")

Dataset 731
time	meso
0.00000

- +

- - + - - + @@ -1093,14 +790,12 @@ kable(pH, caption = "Covariate data")

Dataset 732
time	meso
0.00000

- +

- - + - - + @@ -1124,14 +819,12 @@ kable(pH, caption = "Covariate data")

Dataset 741
time	meso
0.00000

- +

- - + - - + @@ -1157,32 +850,33 @@ kable(pH, caption = "Covariate data")

Dataset 742
time	meso
0.00000

Separate evaluations

Separate evaluations +

In order to obtain suitable starting parameters for the NLHM fits, separate fits of the five models to the data for each soil are generated using the 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", "SFORB", "HS")
-f_sep_const <- mmkin(
-  deg_mods,
-  meso_ds,
-  error_model = "const",
-  cluster = cl,
-  quiet = TRUE)

status(f_sep_const[, 1:5]) |> kable()

- - - +

+deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB", "HS")
+f_sep_const <- mmkin(
+  deg_mods,
+  meso_ds,
+  error_model = "const",
+  cluster = cl,
+  quiet = TRUE)

+status(f_sep_const[, 1:5]) |> kable()

+ - - + @@ -1226,26 +920,26 @@ f_sep_const <- mmkin(

	Richmond	Richmond 2	ERTC	Toulouse	Picket Piece
SFO

status(f_sep_const[, 6:18]) |> kable()

- +

+status(f_sep_const[, 6:18]) |> kable()

--------------++++++++++++++ - - + @@ -1260,8 +954,7 @@ f_sep_const <- mmkin( - - + @@ -1348,19 +1041,19 @@ f_sep_const <- mmkin(

In the tables above, OK indicates convergence and C indicates failure to converge. Most separate fits with constant variance converged, with the exception of two FOMC fits, one SFORB fit and one HS fit.

f_sep_tc <- update(f_sep_const, error_model = "tc")

status(f_sep_tc[, 1:5]) |> kable()

	721	722	732	741	742
SFO

- - +

+f_sep_tc <- update(f_sep_const, error_model = "tc")

+status(f_sep_tc[, 1:5]) |> kable()

+ - - + @@ -1404,26 +1097,26 @@ the exception of two FOMC fits, one SFORB fit and one HS fit.

	Richmond	Richmond 2	ERTC	Toulouse	Picket Piece
SFO

status(f_sep_tc[, 6:18]) |> kable()

- +

+status(f_sep_tc[, 6:18]) |> kable()

--------------++++++++++++++ - - + @@ -1438,8 +1131,7 @@ the exception of two FOMC fits, one SFORB fit and one HS fit.

- - + @@ -1527,21 +1219,21 @@ the exception of two FOMC fits, one SFORB fit and one HS fit.

converge is larger, with convergence problems appearing for a number of non-SFO fits.

Hierarchical model fits without covariate effect

Hierarchical model fits without covariate effect +

The following code fits hierarchical kinetic models for the ten combinations of the five different degradation models with the two different error models in parallel.

f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cluster = cl)
-status(f_saem_1) |> kable()

	721	722	732	741	742
SFO

- - +

+f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cluster = cl)
+status(f_saem_1) |> kable()

+ - - + @@ -1571,17 +1263,16 @@ status(f_saem_1) |> kable()

	const	tc
SFO

All fits terminate without errors (status OK).

anova(f_saem_1) |> kable(digits = 1)

- - - +

+anova(f_saem_1) |> kable(digits = 1)

+ - - + @@ -1660,20 +1351,19 @@ consistently preferable to the corresponding fits with two-component error for these data. This is confirmed by the fact that the parameter b.1 (the relative standard deviation in the fits obtained with the saemix package), is ill-defined in all fits.

illparms(f_saem_1) |> kable()

	npar	AIC	BIC	Lik
SFO const

+illparms(f_saem_1) |> kable()

---+++ - - + - - + @@ -1705,16 +1395,15 @@ with the saemix package), is ill-defined in all fits.

For obtaining fits with only well-defined random effects, we update the set of fits, excluding random effects that were ill-defined according to the illparms function.

f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
-status(f_saem_2) |> kable()

	const	tc
SFO

- - +

+f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+status(f_saem_2) |> kable()

+ - - + @@ -1744,15 +1433,14 @@ status(f_saem_2) |> kable()

	const	tc
SFO

The updated fits terminate without errors.

illparms(f_saem_2) |> kable()

- - - +

+illparms(f_saem_2) |> kable()

+ - - + @@ -1783,8 +1471,9 @@ status(f_saem_2) |> kable()

	const	tc
SFO

No ill-defined errors remain in the fits with constant variance.

Hierarchical model fits with covariate effect

Hierarchical model fits with covariate effect +

In the following sections, hierarchical fits including a model for the influence of pH on selected degradation parameters are shown for all parent models. Constant variance is selected as the error model based on @@ -1793,20 +1482,21 @@ in the fits without pH influence are excluded. A potential influence of the soil pH is only included for parameters with a well-defined random effect, because experience has shown that only for such parameters a significant pH effect could be found.

SFO

sfo_pH <- saem(f_sep_const["SFO", ], no_random_effect = "meso_0", covariates = pH,
-  covariate_models = list(log_k_meso ~ pH))

summary(sfo_pH)$confint_trans |> kable(digits = 2)

- - - +

SFO +

+sfo_pH <- saem(f_sep_const["SFO", ], no_random_effect = "meso_0", covariates = pH,
+  covariate_models = list(log_k_meso ~ pH))

+summary(sfo_pH)$confint_trans |> kable(digits = 2)

+ - - + @@ -1844,23 +1534,26 @@ significant pH effect could be found.

beta_pH(log_k_meso). Its confidence interval does not include zero, indicating that the influence of soil pH on the log of the degradation rate constant is significantly greater than zero.

anova(f_saem_2[["SFO", "const"]], sfo_pH, test = TRUE)

+anova(f_saem_2[["SFO", "const"]], sfo_pH, test = TRUE)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
                            npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
-f_saem_2[["SFO", "const"]]    4 797.56 801.12 -394.78                         
+f_saem_2[["SFO", "const"]]    4 797.56 801.12 -394.78                         
 sfo_pH                        5 783.09 787.54 -386.54 16.473  1  4.934e-05 ***
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The comparison with the SFO fit without covariate effect confirms that considering the soil pH improves the model, both by comparison of AIC and BIC and by the likelihood ratio test.

plot(sfo_pH)

+plot(sfo_pH)

Endpoints for a model with covariates are by default calculated for the median of the covariate values. This quantile can be adapted, or a specific covariate value can be given as shown below.

endpoints(sfo_pH)

+endpoints(sfo_pH)

$covariates
       pH
 50% 5.75
@@ -1868,7 +1561,8 @@ specific covariate value can be given as shown below.
 $distimes
          DT50     DT90
 meso 18.52069 61.52441

endpoints(sfo_pH, covariate_quantile = 0.9)

+endpoints(sfo_pH, covariate_quantile = 0.9)

$covariates
       pH
 90% 7.13
@@ -1876,7 +1570,8 @@ meso 18.52069 61.52441

$distimes DT50 DT90 meso 8.237019 27.36278 -

endpoints(sfo_pH, covariates = c(pH = 7.0))

+endpoints(sfo_pH, covariates = c(pH = 7.0))

$covariates
      pH
 User  7
@@ -1885,20 +1580,21 @@ $distimes
         DT50    DT90
 meso 8.89035 29.5331

FOMC

fomc_pH <- saem(f_sep_const["FOMC", ], no_random_effect = "meso_0", covariates = pH,
-  covariate_models = list(log_alpha ~ pH))

summary(fomc_pH)$confint_trans |> kable(digits = 2)

	est.	lower	upper
meso_0

- - +

FOMC +

+fomc_pH <- saem(f_sep_const["FOMC", ], no_random_effect = "meso_0", covariates = pH,
+  covariate_models = list(log_alpha ~ pH))

+summary(fomc_pH)$confint_trans |> kable(digits = 2)

+ - - + @@ -1951,34 +1647,36 @@ that the model with covariate influence is preferable. However, the random effect for alpha is not well-defined any more after inclusion of the covariate effect (the confidence interval of SD.log_alpha includes zero).

illparms(fomc_pH)

[1] "sd(log_alpha)"

+illparms(fomc_pH)

[1] "sd(log_alpha)"

Therefore, the model is updated without this random effect, and no ill-defined parameters remain.

fomc_pH_2 <- update(fomc_pH, no_random_effect = c("meso_0", "log_alpha"))
-illparms(fomc_pH_2)

anova(f_saem_2[["FOMC", "const"]], fomc_pH, fomc_pH_2, test = TRUE)

+fomc_pH_2 <- update(fomc_pH, no_random_effect = c("meso_0", "log_alpha"))
+illparms(fomc_pH_2)

+anova(f_saem_2[["FOMC", "const"]], fomc_pH, fomc_pH_2, test = TRUE)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
                             npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
-f_saem_2[["FOMC", "const"]]    5 783.25 787.71 -386.63                         
+f_saem_2[["FOMC", "const"]]    5 783.25 787.71 -386.63                         
 fomc_pH_2                      6 767.49 772.83 -377.75 17.762  1  2.503e-05 ***
 fomc_pH                        7 770.07 776.30 -378.04  0.000  1          1    
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model comparison indicates that including pH dependence significantly improves the fit, and that the reduced model with covariate influence results in the most preferable FOMC fit.

summary(fomc_pH_2)$confint_trans |> kable(digits = 2)

	est.	lower	upper
meso_0

- - +

+summary(fomc_pH_2)$confint_trans |> kable(digits = 2)

+ - - + @@ -2018,9 +1716,11 @@ results in the most preferable FOMC fit.

	est.	lower	upper
meso_0

plot(fomc_pH_2)

endpoints(fomc_pH_2)

+plot(fomc_pH_2)

+endpoints(fomc_pH_2)

$covariates
       pH
 50% 5.75
@@ -2028,7 +1728,8 @@ results in the most preferable FOMC fit.
 $distimes
          DT50     DT90 DT50back
 meso 17.30248 82.91343 24.95943

endpoints(fomc_pH_2, covariates = c(pH = 7))

+endpoints(fomc_pH_2, covariates = c(pH = 7))

$covariates
      pH
 User  7
@@ -2037,21 +1738,21 @@ $distimes
          DT50     DT90 DT50back
 meso 6.986239 27.02927 8.136621

DFOP

DFOP +

In the DFOP fits without covariate effects, random effects for two degradation parameters (k2 and g) were identifiable.

summary(f_saem_2[["DFOP", "const"]])$confint_trans |> kable(digits = 2)

- - - +

+summary(f_saem_2[["DFOP", "const"]])$confint_trans |> kable(digits = 2)

+ - - + @@ -2101,23 +1802,23 @@ identifiable.

excluding the same random effects as in the refined DFOP fit without covariate influence, and including covariate models for the two identifiable parameters k2 and g.

dfop_pH <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
-  covariates = pH,
-  covariate_models = list(log_k2 ~ pH, g_qlogis ~ pH))

+dfop_pH <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
+  covariates = pH,
+  covariate_models = list(log_k2 ~ pH, g_qlogis ~ pH))

The corresponding parameters for the influence of soil pH are beta_pH(log_k2) for the influence of soil pH on k2, and beta_pH(g_qlogis) for its influence on g.

summary(dfop_pH)$confint_trans |> kable(digits = 2)

	est.	lower	upper
meso_0

- - +

+summary(dfop_pH)$confint_trans |> kable(digits = 2)

+ - - + @@ -2175,49 +1876,55 @@ identifiable parameters k2 and g.

	est.	lower	upper
meso_0

illparms(dfop_pH)

[1] "sd(g_qlogis)"

+illparms(dfop_pH)

[1] "sd(g_qlogis)"

Confidence intervals for neither of them include zero, indicating a significant difference from zero. However, the random effect for g is now ill-defined. The fit is updated without this ill-defined random effect.

dfop_pH_2 <- update(dfop_pH,
-  no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
-illparms(dfop_pH_2)

[1] "beta_pH(g_qlogis)"

+dfop_pH_2 <- update(dfop_pH,
+  no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
+illparms(dfop_pH_2)

[1] "beta_pH(g_qlogis)"

Now, the slope parameter for the pH effect on g is ill-defined. Therefore, another attempt is made without the corresponding covariate model.

dfop_pH_3 <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
-  covariates = pH,
-  covariate_models = list(log_k2 ~ pH))
-illparms(dfop_pH_3)

[1] "sd(g_qlogis)"

+dfop_pH_3 <- saem(f_sep_const["DFOP", ], no_random_effect = c("meso_0", "log_k1"),
+  covariates = pH,
+  covariate_models = list(log_k2 ~ pH))
+illparms(dfop_pH_3)

[1] "sd(g_qlogis)"

As the random effect for g is again ill-defined, the fit is repeated without it.

dfop_pH_4 <- update(dfop_pH_3, no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
-illparms(dfop_pH_4)

+dfop_pH_4 <- update(dfop_pH_3, no_random_effect = c("meso_0", "log_k1", "g_qlogis"))
+illparms(dfop_pH_4)

While no ill-defined parameters remain, model comparison suggests that the previous model dfop_pH_2 with two pH dependent parameters is preferable, based on information criteria as well as based on the likelihood ratio test.

anova(f_saem_2[["DFOP", "const"]], dfop_pH, dfop_pH_2, dfop_pH_3, dfop_pH_4)

+anova(f_saem_2[["DFOP", "const"]], dfop_pH, dfop_pH_2, dfop_pH_3, dfop_pH_4)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
                             npar    AIC    BIC     Lik
-f_saem_2[["DFOP", "const"]]    7 782.94 789.18 -384.47
+f_saem_2[["DFOP", "const"]]    7 782.94 789.18 -384.47
 dfop_pH_4                      7 767.35 773.58 -376.68
 dfop_pH_2                      8 765.14 772.26 -374.57
 dfop_pH_3                      8 769.00 776.12 -376.50
 dfop_pH                        9 769.10 777.11 -375.55

anova(dfop_pH_2, dfop_pH_4, test = TRUE)

+anova(dfop_pH_2, dfop_pH_4, test = TRUE)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
           npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)  
 dfop_pH_4    7 767.35 773.58 -376.68                       
 dfop_pH_2    8 765.14 772.26 -374.57 4.2153  1    0.04006 *
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

When focussing on parameter identifiability using the test if the confidence interval includes zero, dfop_pH_4 would still be the preferred model. However, it should be kept in mind that parameter @@ -2226,9 +1933,11 @@ likelihood. As the confidence interval of the random effect for g only marginally includes zero, it is suggested that this is acceptable, and that dfop_pH_2 can be considered the most preferable model.

plot(dfop_pH_2)

endpoints(dfop_pH_2)

+plot(dfop_pH_2)

+endpoints(dfop_pH_2)

$covariates
       pH
 50% 5.75
@@ -2236,7 +1945,8 @@ most preferable model.
 $distimes
          DT50     DT90 DT50back  DT50_k1  DT50_k2
 meso 18.36876 73.51841 22.13125 4.191901 23.98672

endpoints(dfop_pH_2, covariates = c(pH = 7))

+endpoints(dfop_pH_2, covariates = c(pH = 7))

$covariates
      pH
 User  7
@@ -2245,21 +1955,22 @@ $distimes
          DT50     DT90 DT50back  DT50_k1  DT50_k2
 meso 8.346428 28.34437 8.532507 4.191901 8.753618

SFORB

sforb_pH <- saem(f_sep_const["SFORB", ], no_random_effect = c("meso_free_0", "log_k_meso_free_bound"),
-  covariates = pH,
-  covariate_models = list(log_k_meso_free ~ pH, log_k_meso_bound_free ~ pH))

summary(sforb_pH)$confint_trans |> kable(digits = 2)

- - - +

SFORB +

+sforb_pH <- saem(f_sep_const["SFORB", ], no_random_effect = c("meso_free_0", "log_k_meso_free_bound"),
+  covariates = pH,
+  covariate_models = list(log_k_meso_free ~ pH, log_k_meso_bound_free ~ pH))

+summary(sforb_pH)$confint_trans |> kable(digits = 2)

+ - - + @@ -2325,39 +2036,41 @@ effect on this parameter (SD.log_k_meso_bound_free) includes zero.

Using the illparms function, these ill-defined parameters can be found more conveniently.

illparms(sforb_pH)

[1] "sd(log_k_meso_bound_free)"      "beta_pH(log_k_meso_bound_free)"

+illparms(sforb_pH)

[1] "sd(log_k_meso_bound_free)"      "beta_pH(log_k_meso_bound_free)"

To remove the ill-defined parameters, a second variant of the SFORB model with pH influence is fitted. No ill-defined parameters remain.

sforb_pH_2 <- update(sforb_pH,
-  no_random_effect = c("meso_free_0", "log_k_meso_free_bound", "log_k_meso_bound_free"),
-  covariate_models = list(log_k_meso_free ~ pH))
-illparms(sforb_pH_2)

+sforb_pH_2 <- update(sforb_pH,
+  no_random_effect = c("meso_free_0", "log_k_meso_free_bound", "log_k_meso_bound_free"),
+  covariate_models = list(log_k_meso_free ~ pH))
+illparms(sforb_pH_2)

The model comparison of the SFORB fits includes the refined model without covariate effect, and both versions of the SFORB fit with covariate effect.

anova(f_saem_2[["SFORB", "const"]], sforb_pH, sforb_pH_2, test = TRUE)

+anova(f_saem_2[["SFORB", "const"]], sforb_pH, sforb_pH_2, test = TRUE)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
                              npar    AIC    BIC     Lik   Chisq Df Pr(>Chisq)  
-f_saem_2[["SFORB", "const"]]    7 783.40 789.63 -384.70                        
+f_saem_2[["SFORB", "const"]]    7 783.40 789.63 -384.70                        
 sforb_pH_2                      7 770.94 777.17 -378.47 12.4616  0             
 sforb_pH                        9 768.81 776.83 -375.41  6.1258  2    0.04675 *
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The first model including pH influence is preferable based on information criteria and the likelihood ratio test. However, as it is not fully identifiable, the second model is selected.

summary(sforb_pH_2)$confint_trans |> kable(digits = 2)

	est.	lower	upper
meso_free_0

- - +

+summary(sforb_pH_2)$confint_trans |> kable(digits = 2)

+ - - + @@ -2403,9 +2116,11 @@ not fully identifiable, the second model is selected.

	est.	lower	upper
meso_free_0

plot(sforb_pH_2)

endpoints(sforb_pH_2)

+plot(sforb_pH_2)

+endpoints(sforb_pH_2)

$covariates
       pH
 50% 5.75
@@ -2421,7 +2136,8 @@ $SFORB
 $distimes
          DT50     DT90 DT50back DT50_meso_b1 DT50_meso_b2
 meso 16.86549 73.15824 22.02282     7.119554     26.33839

endpoints(sforb_pH_2, covariates = c(pH = 7))

+endpoints(sforb_pH_2, covariates = c(pH = 7))

$covariates
      pH
 User  7
@@ -2438,21 +2154,22 @@ $distimes
          DT50     DT90 DT50back DT50_meso_b1 DT50_meso_b2
 meso 7.932495 36.93311 11.11797     5.205671        18.26

HS

hs_pH <- saem(f_sep_const["HS", ], no_random_effect = c("meso_0"),
-  covariates = pH,
-  covariate_models = list(log_k1 ~ pH, log_k2 ~ pH, log_tb ~ pH))

summary(hs_pH)$confint_trans |> kable(digits = 2)

- - - +

HS +

+hs_pH <- saem(f_sep_const["HS", ], no_random_effect = c("meso_0"),
+  covariates = pH,
+  covariate_models = list(log_k1 ~ pH, log_k2 ~ pH, log_tb ~ pH))

+summary(hs_pH)$confint_trans |> kable(digits = 2)

+ - - + @@ -2522,37 +2239,39 @@ meso 7.932495 36.93311 11.11797 5.205671 18.26

	est.	lower	upper
meso_0

illparms(hs_pH)

[1] "sd(log_tb)"      "beta_pH(log_tb)"

+illparms(hs_pH)

[1] "sd(log_tb)"      "beta_pH(log_tb)"

According to the output of the illparms function, the random effect on the break time tb cannot reliably be quantified, neither can the influence of soil pH on tb. The fit is repeated without the corresponding covariate model, and no ill-defined parameters remain.

hs_pH_2 <- update(hs_pH, covariate_models = list(log_k1 ~ pH, log_k2 ~ pH))
-illparms(hs_pH_2)

+hs_pH_2 <- update(hs_pH, covariate_models = list(log_k1 ~ pH, log_k2 ~ pH))
+illparms(hs_pH_2)

Model comparison confirms that this model is preferable to the fit without covariate influence, and also to the first version with covariate influence.

anova(f_saem_2[["HS", "const"]], hs_pH, hs_pH_2, test = TRUE)

+anova(f_saem_2[["HS", "const"]], hs_pH, hs_pH_2, test = TRUE)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
                           npar    AIC    BIC     Lik  Chisq Df Pr(>Chisq)    
-f_saem_2[["HS", "const"]]    8 780.08 787.20 -382.04                         
+f_saem_2[["HS", "const"]]    8 780.08 787.20 -382.04                         
 hs_pH_2                     10 766.47 775.37 -373.23 17.606  2  0.0001503 ***
 hs_pH                       11 769.80 779.59 -373.90  0.000  1  1.0000000    
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(hs_pH_2)$confint_trans |> kable(digits = 2)

- - - +Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +

+summary(hs_pH_2)$confint_trans |> kable(digits = 2)

+ - - + @@ -2616,9 +2335,11 @@ Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.&

	est.	lower	upper
meso_0

plot(hs_pH_2)

endpoints(hs_pH_2)

+plot(hs_pH_2)

+endpoints(hs_pH_2)

$covariates
       pH
 50% 5.75
@@ -2626,7 +2347,8 @@ Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.&
 $distimes
          DT50     DT90 DT50back  DT50_k1  DT50_k2
 meso 14.68725 82.45287 24.82079 14.68725 29.29299

endpoints(hs_pH_2, covariates = c(pH = 7))

+endpoints(hs_pH_2, covariates = c(pH = 7))

$covariates
      pH
 User  7
@@ -2635,11 +2357,13 @@ $distimes
          DT50     DT90 DT50back  DT50_k1  DT50_k2
 meso 8.298536 38.85371 11.69613 8.298536 15.71561

Comparison across parent models

Comparison across parent models +

After model reduction for all models with pH influence, they are compared with each other.

anova(sfo_pH, fomc_pH_2, dfop_pH_2, dfop_pH_4, sforb_pH_2, hs_pH_2)

+anova(sfo_pH, fomc_pH_2, dfop_pH_2, dfop_pH_4, sforb_pH_2, hs_pH_2)

Data: 116 observations of 1 variable(s) grouped in 18 datasets
 
            npar    AIC    BIC     Lik
@@ -2654,7 +2378,8 @@ hs_pH_2      10 766.47 775.37 -373.23

selected as the best fit.

The endpoints resulting from this model are listed below. Please refer to the Appendix for a detailed listing.

endpoints(dfop_pH_2)

+endpoints(dfop_pH_2)

$covariates
       pH
 50% 5.75
@@ -2662,7 +2387,8 @@ refer to the Appendix for a detailed listing.
 $distimes
          DT50     DT90 DT50back  DT50_k1  DT50_k2
 meso 18.36876 73.51841 22.13125 4.191901 23.98672

endpoints(dfop_pH_2, covariates = c(pH = 7))

+endpoints(dfop_pH_2, covariates = c(pH = 7))

$covariates
      pH
 User  7
@@ -2672,34 +2398,40 @@ $distimes
 meso 8.346428 28.34437 8.532507 4.191901 8.753618

Conclusions

Conclusions +

These evaluations demonstrate that covariate effects can be included for all types of parent degradation models. These models can then be further refined to make them fully identifiable.

Appendix

Hierarchical fit listings

Fits without covariate effects

Appendix +

Hierarchical fit listings +

Fits without covariate effects +

Fits with covariate effects

Fits with covariate effects +

Session info

R version 4.4.2 (2024-10-31)
+
+Session info
+
+R version 4.5.0 (2025-04-11)
 Platform: x86_64-pc-linux-gnu
 Running under: Debian GNU/Linux 12 (bookworm)
 
 Matrix products: default
 BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
-LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0
+LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0  LAPACK version 3.11.0
 
 locale:
  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
@@ -2717,74 +2449,51 @@ attached base packages:
 [8] base     
 
 other attached packages:
-[1] saemix_3.3      npde_3.5        knitr_1.49      mkin_1.2.9     
-[5] rmarkdown_2.29  nvimcom_0.9-167
+[1] rmarkdown_2.29  nvimcom_0.9-167 saemix_3.3      npde_3.5       
+[5] knitr_1.49      mkin_1.2.10    
 
 loaded via a namespace (and not attached):
- [1] gtable_0.3.6      jsonlite_1.8.9    dplyr_1.1.4       compiler_4.4.2   
- [5] tidyselect_1.2.1  colorout_1.3-2    gridExtra_2.3     jquerylib_0.1.4  
- [9] scales_1.3.0      readxl_1.4.3      yaml_2.3.10       fastmap_1.2.0    
-[13] lattice_0.22-6    ggplot2_3.5.1     R6_2.5.1          generics_0.1.3   
-[17] lmtest_0.9-40     MASS_7.3-61       tibble_3.2.1      munsell_0.5.1    
-[21] bslib_0.8.0       pillar_1.9.0      rlang_1.1.4       utf8_1.2.4       
-[25] cachem_1.1.0      xfun_0.49         sass_0.4.9        cli_3.6.3        
-[29] magrittr_2.0.3    digest_0.6.37     grid_4.4.2        mclust_6.1.1     
-[33] lifecycle_1.0.4   nlme_3.1-166      vctrs_0.6.5       evaluate_1.0.1   
-[37] glue_1.8.0        cellranger_1.1.0  codetools_0.2-20  zoo_1.8-12       
-[41] fansi_1.0.6       colorspace_2.1-1  tools_4.4.2       pkgconfig_2.0.3  
-[45] htmltools_0.5.8.1
+ [1] gtable_0.3.6      jsonlite_1.9.0    dplyr_1.1.4       compiler_4.5.0   
+ [5] tidyselect_1.2.1  gridExtra_2.3     jquerylib_0.1.4   systemfonts_1.2.1
+ [9] scales_1.3.0      textshaping_1.0.0 readxl_1.4.4      yaml_2.3.10      
+[13] fastmap_1.2.0     lattice_0.22-6    ggplot2_3.5.1     R6_2.6.1         
+[17] generics_0.1.3    lmtest_0.9-40     MASS_7.3-65       htmlwidgets_1.6.4
+[21] tibble_3.2.1      desc_1.4.3        munsell_0.5.1     bslib_0.9.0      
+[25] pillar_1.10.1     rlang_1.1.5       cachem_1.1.0      xfun_0.51        
+[29] fs_1.6.5          sass_0.4.9        cli_3.6.4         pkgdown_2.1.1    
+[33] magrittr_2.0.3    digest_0.6.37     grid_4.5.0        mclust_6.1.1     
+[37] lifecycle_1.0.4   nlme_3.1-168      vctrs_0.6.5       evaluate_1.0.3   
+[41] glue_1.8.0        cellranger_1.1.0  codetools_0.2-20  ragg_1.3.3       
+[45] zoo_1.8-13        colorspace_2.1-1  tools_4.5.0       pkgconfig_2.0.3  
+[49] htmltools_0.5.8.1

Hardware info

Hardware info +

CPU model: AMD Ryzen 9 7950X 16-Core Processor

MemTotal:       64927788 kB

MemTotal:       64927780 kB

- - - - +

Testing covariate modelling in hierarchical -parent degradation kinetics with residue data on mesotrione

Johannes Ranke

Last change on 4 August 2023, last compiled on 13 -Februar 2025

Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione

Johannes +Ranke

Last change on 4 August 2023, +last compiled on 12 Mai 2025

Introduction

Introduction +

Test data

Test data +

Separate evaluations

Separate evaluations +

Hierarchical model fits without covariate effect

Hierarchical model fits without covariate effect +

Hierarchical model fits with covariate effect

Hierarchical model fits with covariate effect +

SFO

SFO +

FOMC

FOMC +

DFOP

DFOP +

SFORB

SFORB +

HS

HS +

Comparison across parent models

Comparison across parent models +

Conclusions

Conclusions +

Appendix

Hierarchical fit listings

Fits without covariate effects

Appendix +

Hierarchical fit listings +

Fits without covariate effects +

Fits with covariate effects

Fits with covariate effects +

Session info

Session info +

Hardware info

Hardware info +