---
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()
```