From 9abab1e2d4385039b01ad3dc0d9c5966bbe94fee Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Aug 2023 17:59:42 +0200 Subject: Update static docs --- docs/articles/prebuilt/2023_mesotrione_parent.html | 2560 ++++++++++++++++++++ .../figure-html/unnamed-chunk-14-1.png | Bin 0 -> 198820 bytes .../figure-html/unnamed-chunk-19-1.png | Bin 0 -> 195998 bytes .../figure-html/unnamed-chunk-25-1.png | Bin 0 -> 199089 bytes .../figure-html/unnamed-chunk-30-1.png | Bin 0 -> 196801 bytes .../figure-html/unnamed-chunk-8-1.png | Bin 0 -> 195237 bytes 6 files changed, 2560 insertions(+) create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent.html create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png create mode 100644 docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png (limited to 'docs/articles/prebuilt') diff --git a/docs/articles/prebuilt/2023_mesotrione_parent.html b/docs/articles/prebuilt/2023_mesotrione_parent.html new file mode 100644 index 00000000..b233fc3c --- /dev/null +++ b/docs/articles/prebuilt/2023_mesotrione_parent.html @@ -0,0 +1,2560 @@ + + + + + + + +Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione • mkin + + + + + + + + + + + + +
+
+ + + + +
+
+
+

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

+

Johannes +Ranke

+ +

Last change on 4 August 2023, +last compiled on 9 August 2023

+ + Source: vignettes/prebuilt/2023_mesotrione_parent.rmd + + +
+ + + +
+

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 +the influence of covariates like soil pH on different degradation +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.5, 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.

+

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.

+
+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")
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Covariate data
pH
Richmond6.2
Richmond 26.2
ERTC6.4
Toulouse7.7
Picket Piece7.1
7215.6
7225.7
7235.4
7244.8
7255.8
7275.1
7285.9
7295.6
7305.3
7316.1
7325.0
7415.7
7427.2
+
+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))
+}
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Richmond
timemeso
0.00000091.00
1.17905086.70
3.53714973.60
7.07429961.50
10.61144855.70
15.32764747.70
17.68574739.50
24.76004629.80
35.37149419.60
68.3848895.67
0.00000097.90
1.17905096.40
3.53714989.10
7.07429974.40
10.61144857.40
15.32764746.30
18.86479735.50
27.11814627.20
35.37149419.10
74.2801386.50
108.4725823.40
142.6650272.20
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Richmond 2
timemeso
0.00000096.0
2.42200482.4
5.65134371.2
8.07334853.1
11.30268748.5
16.95403033.4
22.60537324.2
45.21074611.9
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset ERTC
timemeso
0.00000099.9
2.75519380.0
6.42878242.1
9.18397550.1
12.85756528.4
19.28634739.8
25.71513029.9
51.4302592.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Toulouse
timemeso
0.00000096.8
2.89798363.3
6.76196022.3
9.65994216.6
13.52391916.1
20.28587917.2
27.0478381.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset Picket Piece
timemeso
0.000000102.0
2.84119573.7
6.62945435.5
9.47064931.8
13.25890918.0
19.8883643.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 721
timemeso
0.0000086.4
11.2436661.4
22.4873349.8
33.7309941.0
44.9746635.1
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 722
timemeso
0.0000090.3
11.2436652.1
22.4873337.4
33.7309921.2
44.9746614.3
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 723
timemeso
0.0000089.3
11.2436670.8
22.4873351.1
33.7309942.7
44.9746626.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 724
timemeso
0.00000089.4
9.00820865.2
18.01641555.8
27.02462346.0
36.03283141.7
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 725
timemeso
0.0000089.0
10.9905835.4
21.9811618.6
32.9717411.6
43.962327.6
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 727
timemeso
0.0000091.3
10.9610463.2
21.9220951.1
32.8831342.0
43.8441740.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 728
timemeso
0.0000091.8
11.2436643.6
22.4873322.0
33.7309915.9
44.974668.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 729
timemeso
0.0000091.6
11.2436660.5
22.4873343.5
33.7309928.4
44.9746620.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 730
timemeso
0.0000092.7
11.0744658.9
22.1489344.0
33.2233946.0
44.2978529.3
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 731
timemeso
0.0000092.1
11.2436664.4
22.4873345.3
33.7309933.6
44.9746623.5
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 732
timemeso
0.0000090.3
11.2436658.2
22.4873340.1
33.7309933.1
44.9746625.8
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 741
timemeso
0.0000090.3
10.8471268.7
21.6942458.0
32.5413652.2
43.3884848.0
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dataset 742
timemeso
0.0000092.0
11.2436660.9
22.4873336.2
33.7309918.3
44.974668.7
+
+
+
+

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()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
RichmondRichmond 2ERTCToulousePicket Piece
SFOOKOKOKOKOK
FOMCOKOKOKOKC
DFOPOKOKOKOKOK
SFORBOKOKOKOKOK
HSOKOKCOKOK
+
+status(f_sep_const[, 6:18]) |> kable()
+ ++++++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
721722723724725727728729730731732741742
SFOOKOKOKOKOKOKOKOKOKOKOKOKOK
FOMCOKOKCOKOKOKOKOKOKOKOKOKOK
DFOPOKOKOKOKOKOKOKOKOKOKOKOKOK
SFORBOKOKOKOKOKOKOKCOKOKOKOKOK
HSOKOKOKOKOKOKOKOKOKOKOKOKOK
+

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()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
RichmondRichmond 2ERTCToulousePicket Piece
SFOOKOKOKOKOK
FOMCOKOKOKOKOK
DFOPCOKOKOKOK
SFORBOKOKOKOKOK
HSOKOKCOKOK
+
+status(f_sep_tc[, 6:18]) |> kable()
+ ++++++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
721722723724725727728729730731732741742
SFOOKOKOKOKOKOKOKOKOKOKOKOKOK
FOMCOKOKCOKCCOKCOKCOKCOK
DFOPCOKOKOKCOKOKOKOKCOKCOK
SFORBCOKOKOKCOKOKCOKOKOKCOK
HSOKOKOKOKOKOKOKOKOKCOKOKOK
+

With the two-component error model, the set of fits that did not +converge is larger, with convergence problems appearing for a number of +non-SFO fits.

+
+
+

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()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOOKOK
FOMCOKOK
DFOPOKOK
SFORBOKOK
HSOKOK
+

All fits terminate without errors (status OK).

+
+anova(f_saem_1) |> kable(digits = 1)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
nparAICBICLik
SFO const5800.0804.5-395.0
SFO tc6801.9807.2-394.9
FOMC const7787.4793.6-386.7
FOMC tc8788.9796.1-386.5
DFOP const9787.6795.6-384.8
SFORB const9787.4795.4-384.7
HS const9781.9789.9-382.0
DFOP tc10787.4796.3-383.7
SFORB tc10795.8804.7-387.9
HS tc10783.7792.7-381.9
+

The model comparisons show that the fits with constant variance are +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()
+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOsd(meso_0)sd(meso_0), b.1
FOMCsd(meso_0), sd(log_beta)sd(meso_0), sd(log_beta), b.1
DFOPsd(meso_0), sd(log_k1)sd(meso_0), sd(g_qlogis), b.1
SFORBsd(meso_free_0), sd(log_k_meso_free_bound)sd(meso_free_0), sd(log_k_meso_free_bound), b.1
HSsd(meso_0)sd(meso_0), b.1
+

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()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOOKOK
FOMCOKOK
DFOPOKOK
SFORBOKOK
HSOKOK
+

The updated fits terminate without errors.

+
+illparms(f_saem_2) |> kable()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
consttc
SFOb.1
FOMCb.1
DFOPb.1
SFORBb.1
HS
+

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

+
+
+

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 +the fits without covariate effects. Random effects that were ill-defined +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)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_091.3589.2793.43
log_k_meso-6.66-7.97-5.35
beta_pH(log_k_meso)0.590.370.81
a.15.484.716.24
SD.log_k_meso0.350.230.47
+

The parameter showing the pH influence in the above table is +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)
+
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                         
+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
+

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

+

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)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90
+meso 18.52069 61.52441
+
+endpoints(sfo_pH, covariate_quantile = 0.9)
+
$covariates
+      pH
+90% 7.13
+
+$distimes
+         DT50     DT90
+meso 8.237019 27.36278
+
+endpoints(sfo_pH, covariates = c(pH = 7.0))
+
$covariates
+     pH
+User  7
+
+$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.lowerupper
meso_092.8490.7594.93
log_alpha-2.21-3.49-0.92
beta_pH(log_alpha)0.580.370.79
log_beta4.213.444.99
a.15.034.325.73
SD.log_alpha0.00-23.7723.78
SD.log_beta0.370.010.74
+

As in the case of SFO, the confidence interval of the slope parameter +(here beta_pH(log_alpha)) quantifying the influence of soil +pH does not include zero, and the model comparison clearly indicates +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)"
+

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

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.lowerupper
meso_093.0590.9895.13
log_alpha-2.91-4.18-1.63
beta_pH(log_alpha)0.660.440.87
log_beta3.953.294.62
a.14.984.285.68
SD.log_beta0.400.260.54
+
+plot(fomc_pH_2)
+

+
+endpoints(fomc_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$distimes
+         DT50     DT90 DT50back
+meso 17.30248 82.91343 24.95943
+
+endpoints(fomc_pH_2, covariates = c(pH = 7))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back
+meso 6.986239 27.02927 8.136621
+
+
+

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)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.6191.5895.63
log_k1-1.53-2.27-0.79
log_k2-3.42-3.73-3.11
g_qlogis-1.67-2.57-0.77
a.14.744.025.45
SD.log_k20.600.380.81
SD.g_qlogis0.940.331.54
+

A fit with pH dependent degradation parameters was obtained by +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))
+

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.lowerupper
meso_092.8490.8594.84
log_k1-2.82-3.09-2.54
log_k2-11.48-15.32-7.64
beta_pH(log_k2)1.310.691.92
g_qlogis3.130.475.80
beta_pH(g_qlogis)-0.57-1.04-0.09
a.14.964.265.65
SD.log_k20.760.471.05
SD.g_qlogis0.01-9.969.97
+
+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)"
+

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

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

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

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 +confidence intervals are constructed using a simple linearisation of the +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)
+
$covariates
+      pH
+50% 5.75
+
+$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))
+
$covariates
+     pH
+User  7
+
+$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)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_free_093.4291.3295.52
log_k_meso_free-5.37-6.94-3.81
beta_pH(log_k_meso_free)0.420.180.67
log_k_meso_free_bound-3.49-4.92-2.05
log_k_meso_bound_free-9.98-19.22-0.74
beta_pH(log_k_meso_bound_free)1.23-0.212.67
a.14.904.185.63
SD.log_k_meso_free0.350.230.47
SD.log_k_meso_bound_free0.13-1.952.20
+

The confidence interval of +beta_pH(log_k_meso_bound_free) includes zero, indicating +that the influence of soil pH on k_meso_bound_free cannot +reliably be quantified. Also, the confidence interval for the random +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)"
+

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

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

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.lowerupper
meso_free_093.3291.1695.48
log_k_meso_free-6.15-7.43-4.86
beta_pH(log_k_meso_free)0.540.330.75
log_k_meso_free_bound-3.80-5.20-2.40
log_k_meso_bound_free-2.95-4.26-1.64
a.15.084.385.79
SD.log_k_meso_free0.330.220.45
+
+plot(sforb_pH_2)
+

+
+endpoints(sforb_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$ff
+meso_free 
+        1 
+
+$SFORB
+   meso_b1    meso_b2     meso_g 
+0.09735824 0.02631699 0.31602120 
+
+$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))
+
$covariates
+     pH
+User  7
+
+$ff
+meso_free 
+        1 
+
+$SFORB
+   meso_b1    meso_b2     meso_g 
+0.13315233 0.03795988 0.61186191 
+
+$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)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.3391.4795.19
log_k1-5.81-7.27-4.36
beta_pH(log_k1)0.470.230.72
log_k2-6.80-8.76-4.83
beta_pH(log_k2)0.540.210.87
log_tb3.251.255.25
beta_pH(log_tb)-0.10-0.430.23
a.14.493.785.21
SD.log_k10.370.240.51
SD.log_k20.290.100.48
SD.log_tb0.25-0.070.57
+
+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)
+

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)
+
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                         
+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)
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
est.lowerupper
meso_093.3391.5095.15
log_k1-5.68-7.09-4.27
beta_pH(log_k1)0.460.220.69
log_k2-6.61-8.34-4.88
beta_pH(log_k2)0.500.210.79
log_tb2.702.333.08
a.14.453.745.16
SD.log_k10.360.220.49
SD.log_k20.230.020.43
SD.log_tb0.550.250.85
+
+plot(hs_pH_2)
+

+
+endpoints(hs_pH_2)
+
$covariates
+      pH
+50% 5.75
+
+$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))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 8.298536 38.85371 11.69613 8.298536 15.71561
+
+
+

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)
+
Data: 116 observations of 1 variable(s) grouped in 18 datasets
+
+           npar    AIC    BIC     Lik
+sfo_pH        5 783.09 787.54 -386.54
+fomc_pH_2     6 767.49 772.83 -377.75
+dfop_pH_4     7 767.35 773.58 -376.68
+sforb_pH_2    7 770.94 777.17 -378.47
+dfop_pH_2     8 765.14 772.26 -374.57
+hs_pH_2      10 766.47 775.37 -373.23
+

The DFOP model with pH influence on k2 and +g and a random effect only on k2 is finally +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)
+
$covariates
+      pH
+50% 5.75
+
+$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))
+
$covariates
+     pH
+User  7
+
+$distimes
+         DT50     DT90 DT50back  DT50_k1  DT50_k2
+meso 8.346428 28.34437 8.532507 4.191901 8.753618
+
+
+
+

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 +

+ +
+
+

Fits with covariate effects +

+ +
+
+
+

Session info +

+
R version 4.3.1 (2023-06-16)
+Platform: x86_64-pc-linux-gnu (64-bit)
+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
+
+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       
+
+time zone: Europe/Berlin
+tzcode source: system (glibc)
+
+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.43 mkin_1.2.5
+
+loaded via a namespace (and not attached):
+ [1] sass_0.4.6        utf8_1.2.3        generics_0.1.3    stringi_1.7.12   
+ [5] lattice_0.20-45   digest_0.6.31     magrittr_2.0.3    evaluate_0.21    
+ [9] grid_4.3.1        fastmap_1.1.1     cellranger_1.1.0  rprojroot_2.0.3  
+[13] jsonlite_1.8.5    mclust_6.0.0      gridExtra_2.3     purrr_1.0.1      
+[17] fansi_1.0.4       scales_1.2.1      codetools_0.2-19  textshaping_0.3.6
+[21] jquerylib_0.1.4   cli_3.6.1         rlang_1.1.1       munsell_0.5.0    
+[25] cachem_1.0.8      yaml_2.3.7        tools_4.3.1       memoise_2.0.1    
+[29] dplyr_1.1.2       colorspace_2.1-0  ggplot2_3.4.2     vctrs_0.6.2      
+[33] R6_2.5.1          zoo_1.8-12        lifecycle_1.0.3   stringr_1.5.0    
+[37] fs_1.6.2          ragg_1.2.5        pkgconfig_2.0.3   desc_1.4.2       
+[41] pkgdown_2.0.7     bslib_0.4.2       pillar_1.9.0      gtable_0.3.3     
+[45] glue_1.6.2        systemfonts_1.0.4 highr_0.10        xfun_0.39        
+[49] tibble_3.2.1      lmtest_0.9-40     tidyselect_1.2.0  htmltools_0.5.5  
+[53] nlme_3.1-162      rmarkdown_2.22    compiler_4.3.1    readxl_1.4.2
+
+
+

Hardware info +

+
CPU model: Intel(R) Core(TM) i7-4710MQ CPU @ 2.50GHz
+
MemTotal:       12165632 kB
+
+
+
+ + + +
+ + + + +
+ + + + + + + + diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png new file mode 100644 index 00000000..863a48bd Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png new file mode 100644 index 00000000..256b2b68 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png new file mode 100644 index 00000000..59011020 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-25-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png new file mode 100644 index 00000000..f427bc39 Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-30-1.png differ diff --git a/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png new file mode 100644 index 00000000..7c3b460b Binary files /dev/null and b/docs/articles/prebuilt/2023_mesotrione_parent_files/figure-html/unnamed-chunk-8-1.png differ -- cgit v1.2.1