1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
---
title: "Hierarchical kinetic modelling of degradation data"
author:
date:
output: mkin::hierarchical_kinetics
geometry: margin=2cm
---
```{r packages, cache = FALSE, message = FALSE}
library(mkin)
library(knitr)
library(saemix)
library(parallel)
library(readxl)
```
```{r n_cores, cache = FALSE, echo = FALSE}
n_cores <- detectCores()
if (Sys.info()["sysname"] == "Windows") {
cl <- makePSOCKcluster(n_cores)
} else {
cl <- makeForkCluster(n_cores)
}
```
\clearpage
# Introduction
This report shows hierarchical kinetic modelling for ...
The data were obtained from ...
```{r ds}
data_path <- system.file("testdata", "lambda-cyhalothrin_soil_efsa_2014.xlsx", package = "mkin")
ds <- read_spreadsheet(data_path, valid_datasets = c(1:4, 7:13))
covariates <- attr(ds, "covariates")
```
The covariate data are shown below.
```{r results = "asis", dependson = "ds", echo = FALSE}
kable(covariates, caption = "Covariate data for all datasets")
```
\clearpage
The datasets with the residue time series are shown in the tables below. Please
refer to the spreadsheet for details like data sources, treatment of values
below reporting limits and time step normalisation factors.
```{r results = "asis", dependson = "ds", echo = FALSE}
for (ds_name in names(ds)) {
print(
kable(mkin_long_to_wide(ds[[ds_name]]),
caption = paste("Dataset", ds_name),
booktabs = TRUE, row.names = FALSE))
cat("\n\\clearpage\n")
}
```
# Parent only evaluations
The following code performs separate fits of the candidate degradation models
to all datasets using constant variance and the two-component error model.
```{r parent-sep, dependson = "ds"}
parent_deg_mods <- c("SFO", "FOMC", "DFOP", "SFORB")
errmods <- c(const = "constant variance", tc = "two-component error")
parent_sep_const <- mmkin(
parent_deg_mods, ds,
error_model = "const",
cluster = cl, quiet = TRUE)
parent_sep_tc <- update(parent_sep_const, error_model = "tc")
```
To select the parent model, the corresponding hierarchical fits are performed below.
```{r parent-mhmkin, dependson = "parent-sep"}
parent_mhmkin <- mhmkin(list(parent_sep_const, parent_sep_tc), cluster = cl)
status(parent_mhmkin) |> kable()
```
All fits terminate without errors (status OK). The check for ill-defined
parameters shows that not all random effect parameters can be robustly
quantified.
```{r dependson = "parent_mhmkin"}
illparms(parent_mhmkin) |> kable()
```
Therefore, the fits are updated, excluding random effects that were
ill-defined according to the `illparms` function.
```{r parent-mhmkin-refined}
parent_mhmkin_refined <- update(parent_mhmkin,
no_random_effect = illparms(parent_mhmkin))
status(parent_mhmkin_refined) |> kable()
```
The most suitable model is selected based on the AIC.
```{r, dependson = "parent-mhmkin"}
aic_parent <- AIC(parent_mhmkin_refined)
min_aic <- which(aic_parent == min(aic_parent), arr.ind = TRUE)
best_degmod_parent <- rownames(aic_parent)[min_aic[1]]
best_errmod_parent <- colnames(aic_parent)[min_aic[2]]
anova(parent_mhmkin_refined) |> kable(digits = 1)
```
Based on the AIC, the combination of the `r best_degmod_parent` degradation
model with the error model `r errmods[best_errmod_parent]` is identified to
be most suitable for the degradation of the parent. The check below
confirms that no ill-defined parameters remain for this combined model.
```{r dependson = "parent-mhmkin"}
illparms(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]])
```
The corresponding fit is shown below.
```{r parent-best-full, dependson = "parent-mhmkin"}
plot(parent_mhmkin_refined[[best_degmod_parent, best_errmod_parent]])
```
Detailed listings of the parent fits are shown in the Appendix.
\clearpage
# Appendix
## Summaries of saem fits
### Parent fits
```{r listings-parent, results = "asis", echo = FALSE, dependson = "parent_mhmkin"}
for (deg_mod in parent_deg_mods) {
for (err_mod in c("const", "tc")) {
caption <- paste("Hierarchical", deg_mod, "fit with", errmods[err_mod])
tex_listing(parent_mhmkin[[deg_mod, err_mod]], caption)
}
}
```
### Refined parent fits
```{r listings-parent-refined, results = "asis", echo = FALSE, dependson = "parent_mhmkin_refined"}
for (deg_mod in parent_deg_mods) {
for (err_mod in c("const", "tc")) {
caption <- paste("Refined hierarchical", deg_mod, "fit with", errmods[err_mod])
tex_listing(parent_mhmkin_refined[[deg_mod, err_mod]], caption)
}
}
```
## Session info
```{r, echo = FALSE, cache = FALSE}
parallel::stopCluster(cl)
sessionInfo()
```
|
