From df0cff4b829f1abf62f037591a24a8019174dd0a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 18 Nov 2022 08:37:40 +0100 Subject: Pass error.init to saemix_model, show in parplot Due to an oversight, error.init was not really passed to saemix_model in saem.mmkin. The new initial values were reverted to c(1, 1), in order to avoid changing the test results. Initial values for error model parameters are now shown in parplot.multistart. --- R/parplot.R | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) (limited to 'R/parplot.R') diff --git a/R/parplot.R b/R/parplot.R index 627a4eb8..98579779 100644 --- a/R/parplot.R +++ b/R/parplot.R @@ -4,6 +4,10 @@ #' either by the parameters of the run with the highest likelihood, #' or by their medians as proposed in the paper by Duchesne et al. (2021). #' +#' Starting values of degradation model parameters and error model parameters +#' are shown as green circles. The results obtained in the original run +#' are shown as red circles. +#' #' @param object The [multistart] object #' @param llmin The minimum likelihood of objects to be shown #' @param scale By default, scale parameters using the best available fit. @@ -32,7 +36,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best" orig <- attr(object, "orig") orig_parms <- parms(orig) - start_parms <- orig$mean_dp_start + start_degparms <- orig$mean_dp_start all_parms <- parms(object) if (inherits(object, "multistart.saem.mmkin")) { @@ -49,18 +53,18 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best" par(las = 1) if (orig$transformations == "mkin") { - degparm_names_transformed <- names(start_parms) + degparm_names_transformed <- names(start_degparms) degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) orig_parms[degparm_names_transformed] <- backtransform_odeparms( orig_parms[degparm_names_transformed], orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions) - start_parms <- backtransform_odeparms(start_parms, + start_degparms <- backtransform_odeparms(start_degparms, orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions) - degparm_names <- names(start_parms) + degparm_names <- names(start_degparms) names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) @@ -72,6 +76,13 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best" colnames(selected_parms)[1:length(degparm_names)] <- degparm_names } + start_errparms <- orig$so@model@error.init + names(start_errparms) <- orig$so@model@name.sigma + + start_omegaparms <- orig$so@model@omega.init + + start_parms <- c(start_degparms, start_errparms) + scale <- match.arg(scale) parm_scale <- switch(scale, best = selected_parms[which.best(object[selected]), ], -- cgit v1.2.3 From 478c6d5eec4c84b22b43adcbdf36888b302ead00 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 6 Dec 2022 10:33:24 +0100 Subject: Some parplot improvements llquant argument, improved legend text, tests --- NEWS.md | 2 + R/parplot.R | 16 +- log/test.log | 41 +++-- man/parplot.Rd | 7 +- .../_snaps/multistart/parplot-for-dfop-sfo-fit.svg | 169 ++++++++++----------- .../_snaps/multistart/parplot-for-sfo-fit.svg | 4 +- tests/testthat/test_multistart.R | 2 +- vignettes/FOCUS_D.html | 37 ++--- vignettes/FOCUS_L.html | 83 +++++----- 9 files changed, 183 insertions(+), 178 deletions(-) (limited to 'R/parplot.R') diff --git a/NEWS.md b/NEWS.md index 15c45cc9..c7929fbe 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,8 @@ - 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted. +- 'R/parplot.R': Possibility to select the top 'llquant' fraction of the fits for the parameter plots, and improved legend text. + # mkin 1.2.1 (2022-11-19) - '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests. diff --git a/R/parplot.R b/R/parplot.R index 98579779..63306ac2 100644 --- a/R/parplot.R +++ b/R/parplot.R @@ -10,7 +10,10 @@ #' #' @param object The [multistart] object #' @param llmin The minimum likelihood of objects to be shown -#' @param scale By default, scale parameters using the best available fit. +#' @param llquant Fractional value for selecting only the fits with higher +#' likelihoods. Overrides 'llmin'. +#' @param scale By default, scale parameters using the best +#' available fit. #' If 'median', parameters are scaled using the median parameters from all fits. #' @param main Title of the plot #' @param lpos Positioning of the legend. @@ -28,7 +31,8 @@ parplot <- function(object, ...) { #' @rdname parplot #' @export -parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"), +parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, + scale = c("best", "median"), lpos = "bottomleft", main = "", ...) { oldpar <- par(no.readonly = TRUE) @@ -48,6 +52,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best" stop("parplot is only implemented for multistart.saem.mmkin objects") } ll <- sapply(object, llfunc) + if (!is.na(llquant[1])) { + if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified") + llmin <- quantile(ll, 1 - llquant) + } selected <- which(ll > llmin) selected_parms <- all_parms[selected, ] @@ -110,7 +118,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best" legend(lpos, inset = c(0.05, 0.05), bty = "n", pch = 1, col = 3:1, lty = c(NA, NA, 1), legend = c( - "Starting parameters", - "Original run", + "Original start", + "Original results", "Multistart runs")) } diff --git a/log/test.log b/log/test.log index 5764f209..7614b136 100644 --- a/log/test.log +++ b/log/test.log @@ -1,22 +1,22 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.3s] +✔ | 5 | Analytical solutions for coupled models [3.2s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.1s] -✔ | 1 12 | Dimethenamid data from 2018 [31.9s] +✔ | 1 12 | Dimethenamid data from 2018 [32.0s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ──────────────────────────────────────────────────────────────────────────────── -✔ | 14 | Error model fitting [4.9s] +✔ | 14 | Error model fitting [4.8s] ✔ | 5 | Time step normalisation ✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] ✔ | 1 | Fitting the logistic model [0.2s] -✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [41.4s] +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [42.2s] ✔ | 1 11 | Nonlinear mixed-effects models [13.3s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits @@ -26,32 +26,41 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 10 | Special cases of mkinfit calls [0.6s] ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] -✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 12 | Multistart method for saem.mmkin models [47.5s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.3s] -✔ | 9 | Nonlinear mixed-effects models with nlme [9.5s] -✔ | 15 | Plotting [10.2s] +✔ | 3 | Model predictions with mkinpredict [0.4s] +✖ | 1 11 | Multistart method for saem.mmkin models [46.7s] +──────────────────────────────────────────────────────────────────────────────── +Failure ('test_multistart.R:56'): multistart works for saem.mmkin models +Snapshot of `testcase` to 'multistart/parplot-for-dfop-sfo-fit.svg' has changed +Run `testthat::snapshot_review('multistart/')` to review changes +Backtrace: + 1. vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo) + at test_multistart.R:56:2 + 3. testthat::expect_snapshot_file(...) +──────────────────────────────────────────────────────────────────────────────── +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s] +✔ | 9 | Nonlinear mixed-effects models with nlme [9.4s] +✔ | 15 | Plotting [10.1s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [73.0s] +✔ | 1 36 | saemix parent models [73.6s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem Reason: This still takes almost 2.5 minutes although we do not solve ODEs ──────────────────────────────────────────────────────────────────────────────── -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] ✔ | 11 | Processing of residue series -✔ | 10 | Fitting the SFORB model [3.9s] +✔ | 10 | Fitting the SFORB model [3.7s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary [0.2s] -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s] -✔ | 9 | Hypothesis tests [8.4s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s] +✔ | 9 | Hypothesis tests [8.1s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 261.1 s +Duration: 260.9 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) • Fitting with saemix takes around 10 minutes when using deSolve (1) • This still takes almost 2.5 minutes although we do not solve ODEs (1) -[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ] +[ FAIL 1 | WARN 0 | SKIP 3 | PASS 269 ] diff --git a/man/parplot.Rd b/man/parplot.Rd index ac9e02cf..67bf0cc1 100644 --- a/man/parplot.Rd +++ b/man/parplot.Rd @@ -10,6 +10,7 @@ parplot(object, ...) \method{parplot}{multistart.saem.mmkin}( object, llmin = -Inf, + llquant = NA, scale = c("best", "median"), lpos = "bottomleft", main = "", @@ -23,7 +24,11 @@ parplot(object, ...) \item{llmin}{The minimum likelihood of objects to be shown} -\item{scale}{By default, scale parameters using the best available fit. +\item{llquant}{Fractional value for selecting only the fits with higher +likelihoods. Overrides 'llmin'.} + +\item{scale}{By default, scale parameters using the best +available fit. If 'median', parameters are scaled using the median parameters from all fits.} \item{lpos}{Positioning of the legend.} diff --git a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg index 7017908e..b01dac74 100644 --- a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg +++ b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg @@ -25,109 +25,104 @@ - - - - + + + + - - - - - - - - + + + + + + - - - - - + + + + + - + - - - - - - - - - - - - - - + + + + + + + + + + + + - - - - - + + + + + - - - - - - + + + + + + - - - - - + + + + + - - - - - - - - + + + + + + + - - - - - + + + + + - - - - - - + + + + + + - - - - - - + + + + + + - - - - - - + + + + + + - - - - - + + + + + - + @@ -194,8 +189,8 @@ -Starting parameters -Original run +Original start +Original results Multistart runs diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg index a47a585a..c8d4970b 100644 --- a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg +++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg @@ -103,8 +103,8 @@ -Starting parameters -Original run +Original start +Original results Multistart runs diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R index 91ef71f0..dda0ea23 100644 --- a/tests/testthat/test_multistart.R +++ b/tests/testthat/test_multistart.R @@ -50,7 +50,7 @@ test_that("multistart works for saem.mmkin models", { llhist_dfop_sfo <- function() llhist(saem_dfop_sfo_m_multi) parplot_dfop_sfo <- function() parplot(saem_dfop_sfo_m_multi, - ylim = c(0.5, 2)) + ylim = c(0.5, 2), llquant = 0.5) vdiffr::expect_doppelganger("llhist for dfop sfo fit", llhist_dfop_sfo) vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo) diff --git a/vignettes/FOCUS_D.html b/vignettes/FOCUS_D.html index 9f768b06..b8a63a7b 100644 --- a/vignettes/FOCUS_D.html +++ b/vignettes/FOCUS_D.html @@ -299,8 +299,8 @@ pre code { border-radius: 4px; } -.tabset-dropdown > .nav-tabs > li.active:before { - content: ""; +.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before { + content: "\e259"; font-family: 'Glyphicons Halflings'; display: inline-block; padding: 10px; @@ -308,16 +308,9 @@ pre code { } .tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before { - content: ""; - border: none; -} - -.tabset-dropdown > .nav-tabs.nav-tabs-open:before { - content: ""; + content: "\e258"; font-family: 'Glyphicons Halflings'; - display: inline-block; - padding: 10px; - border-right: 1px solid #ddd; + border: none; } .tabset-dropdown > .nav-tabs > li.active { @@ -364,7 +357,7 @@ pre code {

Example evaluation of FOCUS Example Dataset D

Johannes Ranke

-

Last change 31 January 2019 (rebuilt 2022-07-08)

+

Last change 31 January 2019 (rebuilt 2022-12-06)

@@ -438,10 +431,10 @@ print(FOCUS_2006_D)

A comprehensive report of the results is obtained using the summary method for mkinfit objects.

summary(fit)
-
## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Fri Jul  8 15:44:37 2022 
-## Date of summary: Fri Jul  8 15:44:38 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:42 2022 
+## Date of summary: Tue Dec  6 09:39:42 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
@@ -449,7 +442,7 @@ print(FOCUS_2006_D)

 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 401 model solutions performed in 0.13 s
+## Fitted using 401 model solutions performed in 0.158 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -492,11 +485,11 @@ print(FOCUS_2006_D)
## ## Parameter correlation: ## parent_0 log_k_parent log_k_m1 f_parent_qlogis sigma -## parent_0 1.000e+00 5.174e-01 -1.688e-01 -5.471e-01 -1.174e-06 -## log_k_parent 5.174e-01 1.000e+00 -3.263e-01 -5.426e-01 -8.492e-07 -## log_k_m1 -1.688e-01 -3.263e-01 1.000e+00 7.478e-01 8.220e-07 -## f_parent_qlogis -5.471e-01 -5.426e-01 7.478e-01 1.000e+00 1.307e-06 -## sigma -1.174e-06 -8.492e-07 8.220e-07 1.307e-06 1.000e+00 +## parent_0 1.000e+00 5.174e-01 -1.688e-01 -5.471e-01 -1.172e-06 +## log_k_parent 5.174e-01 1.000e+00 -3.263e-01 -5.426e-01 -8.483e-07 +## log_k_m1 -1.688e-01 -3.263e-01 1.000e+00 7.478e-01 8.205e-07 +## f_parent_qlogis -5.471e-01 -5.426e-01 7.478e-01 1.000e+00 1.305e-06 +## sigma -1.172e-06 -8.483e-07 8.205e-07 1.305e-06 1.000e+00 ## ## Backtransformed parameters: ## Confidence intervals for internally transformed parameters are asymmetric. diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html index da6c11fe..b8c9ba0c 100644 --- a/vignettes/FOCUS_L.html +++ b/vignettes/FOCUS_L.html @@ -1373,8 +1373,8 @@ pre code { border-radius: 4px; } -.tabset-dropdown > .nav-tabs > li.active:before { - content: ""; +.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before { + content: "\e259"; font-family: 'Glyphicons Halflings'; display: inline-block; padding: 10px; @@ -1382,16 +1382,9 @@ pre code { } .tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before { - content: ""; - border: none; -} - -.tabset-dropdown > .nav-tabs.nav-tabs-open:before { - content: ""; + content: "\e258"; font-family: 'Glyphicons Halflings'; - display: inline-block; - padding: 10px; - border-right: 1px solid #ddd; + border: none; } .tabset-dropdown > .nav-tabs > li.active { @@ -1517,7 +1510,7 @@ div.tocify {

Example evaluation of FOCUS Laboratory Data L1 to L3

Johannes Ranke

-

Last change 18 May 2022 (rebuilt 2022-09-14)

+

Last change 18 May 2022 (rebuilt 2022-12-06)

@@ -1536,17 +1529,17 @@ FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1)

Since mkin version 0.9-32 (July 2014), we can use shorthand notation like "SFO" for parent only degradation models. The following two lines fit the model and produce the summary report of the model fit. This covers the numerical analysis given in the FOCUS report.

m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE)
 summary(m.L1.SFO)
-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:35 2022 
-## Date of summary: Wed Sep 14 22:28:35 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:45 2022 
+## Date of summary: Tue Dec  6 09:39:45 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 133 model solutions performed in 0.032 s
+## Fitted using 133 model solutions performed in 0.033 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1637,10 +1630,10 @@ summary(m.L1.SFO)

 
## Warning in sqrt(1/diag(V)): NaNs produced

 
## Warning in cov2cor(ans$covar): diag(.) had 0 or NA entries; non-finite result is
 ## doubtful

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:35 2022 
-## Date of summary: Wed Sep 14 22:28:35 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:45 2022 
+## Date of summary: Tue Dec  6 09:39:45 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
@@ -1742,17 +1735,17 @@ plot(m.L2.FOMC, show_residuals = TRUE,
      main = "FOCUS L2 - FOMC")

 

 
summary(m.L2.FOMC, data = FALSE)

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:35 2022 
-## Date of summary: Wed Sep 14 22:28:35 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:45 2022 
+## Date of summary: Tue Dec  6 09:39:45 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 239 model solutions performed in 0.049 s
+## Fitted using 239 model solutions performed in 0.048 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1820,10 +1813,10 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
      main = "FOCUS L2 - DFOP")

 

 
summary(m.L2.DFOP, data = FALSE)

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:36 2022 
-## Date of summary: Wed Sep 14 22:28:36 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:46 2022 
+## Date of summary: Tue Dec  6 09:39:46 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1832,7 +1825,7 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 581 model solutions performed in 0.135 s
+## Fitted using 581 model solutions performed in 0.131 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1920,10 +1913,10 @@ plot(mm.L3)

 
The objects returned by mmkin are arranged like a matrix, with models as a row index and datasets as a column index.

 
We can extract the summary and plot for e.g. the DFOP fit, using square brackets for indexing which will result in the use of the summary and plot functions working on mkinfit objects.

 
summary(mm.L3[["DFOP", 1]])

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:36 2022 
-## Date of summary: Wed Sep 14 22:28:36 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:46 2022 
+## Date of summary: Tue Dec  6 09:39:46 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1932,7 +1925,7 @@ plot(mm.L3)

 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 376 model solutions performed in 0.081 s
+## Fitted using 376 model solutions performed in 0.078 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -2028,17 +2021,17 @@ plot(mm.L4)

 

 
The χ2 error level of 3.3% as well as the plot suggest that the SFO model fits very well. The error level at which the χ2 test passes is slightly lower for the FOMC model. However, the difference appears negligible.

 
summary(mm.L4[["SFO", 1]], data = FALSE)

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:36 2022 
-## Date of summary: Wed Sep 14 22:28:37 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:47 2022 
+## Date of summary: Tue Dec  6 09:39:47 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 142 model solutions performed in 0.034 s
+## Fitted using 142 model solutions performed in 0.03 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -2092,10 +2085,10 @@ plot(mm.L4)

 ##        DT50 DT90
 ## parent  106  352

 
summary(mm.L4[["FOMC", 1]], data = FALSE)

-
## mkin version used for fitting:    1.1.2 
-## R version used for fitting:       4.2.1 
-## Date of fit:     Wed Sep 14 22:28:37 2022 
-## Date of summary: Wed Sep 14 22:28:37 2022 
+
## mkin version used for fitting:    1.2.2 
+## R version used for fitting:       4.2.2 
+## Date of fit:     Tue Dec  6 09:39:47 2022 
+## Date of summary: Tue Dec  6 09:39:47 2022 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
-- 
cgit v1.2.3


From 904ba9668eb76eaae4960e2188134e8c88da07ee Mon Sep 17 00:00:00 2001
From: Johannes Ranke 
Date: Wed, 7 Dec 2022 16:19:54 +0100
Subject: Fix parplot for the case of failed multistart runs

---
 NEWS.md     | 2 ++
 R/parms.R   | 2 +-
 R/parplot.R | 2 +-
 3 files changed, 4 insertions(+), 2 deletions(-)

(limited to 'R/parplot.R')

diff --git a/NEWS.md b/NEWS.md
index 1931142c..4540b517 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,8 @@
 
 - 'R/illparms.R': Also check if confidence intervals for slope parameters in covariate models include zero. Only implemented for fits obtained with the saemix backend.
 
+- 'R/parplot.R': Make the function work also in the case that some of the multistart runs failed.
+
 # mkin 1.2.1 (2022-11-19)
 
 - '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/parms.R b/R/parms.R
index bd4e479b..bb04a570 100644
--- a/R/parms.R
+++ b/R/parms.R
@@ -77,6 +77,6 @@ parms.multistart <- function(object, exclude_failed = TRUE, ...) {
   successful <- which(!is.na(res[, 1]))
   first_success <- successful[1]
   colnames(res) <- names(parms(object[[first_success]]))
-  if (exclude_failed) res <- res[successful, ]
+  if (exclude_failed[1]) res <- res[successful, ]
   return(res)
 }
diff --git a/R/parplot.R b/R/parplot.R
index 63306ac2..e9c18947 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -41,7 +41,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
   orig <- attr(object, "orig")
   orig_parms <- parms(orig)
   start_degparms <- orig$mean_dp_start
-  all_parms <- parms(object)
+  all_parms <- parms(object, exclude_failed = FALSE)
 
   if (inherits(object, "multistart.saem.mmkin")) {
     llfunc <- function(object) {
-- 
cgit v1.2.3


From 886c9ef013124aa954d960c655b349b5340ff154 Mon Sep 17 00:00:00 2001
From: Johannes Ranke 
Date: Mon, 19 Dec 2022 12:31:56 +0100
Subject: Rename template folder, create format

Instead of rmarkdown::pdf_document, mkin::hierarchical_kinetics is used
as a document format in the template. In this way, the template file can
be freed from some R code and yaml options that the average user does
not have to be aware of.
---
 .Rbuildignore                                      |   6 +-
 .gitignore                                         |   6 +-
 GNUmakefile                                        |   4 +-
 NAMESPACE                                          |   2 +
 NEWS.md                                            |   4 +
 R/hierarchical_kinetics.R                          |  39 +++++
 R/parplot.R                                        |   2 +-
 .../rmarkdown/templates/hier/skeleton/skeleton.Rmd | 186 ---------------------
 inst/rmarkdown/templates/hier/template.yaml        |   3 -
 .../hierarchical_kinetics/skeleton/header.tex      |   1 +
 .../hierarchical_kinetics/skeleton/skeleton.Rmd    | 162 ++++++++++++++++++
 .../templates/hierarchical_kinetics/template.yaml  |   3 +
 log/check.log                                      |  21 ++-
 man/hierarchical_kinetics.Rd                       |  29 ++++
 14 files changed, 265 insertions(+), 203 deletions(-)
 create mode 100644 R/hierarchical_kinetics.R
 delete mode 100644 inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd
 delete mode 100644 inst/rmarkdown/templates/hier/template.yaml
 create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex
 create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
 create mode 100644 inst/rmarkdown/templates/hierarchical_kinetics/template.yaml
 create mode 100644 man/hierarchical_kinetics.Rd

(limited to 'R/parplot.R')

diff --git a/.Rbuildignore b/.Rbuildignore
index 0dcec29a..1eb33577 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,9 +9,9 @@
 ^test.R$
 ^mkin_.*\.tar\.gz
 ^mkin.tar$
-^inst/rmarkdown/templates/hier/skeleton/skeleton.pdf$
-^inst/rmarkdown/templates/hier/skeleton/skeleton_cache
-^inst/rmarkdown/templates/hier/skeleton/skeleton_files
+^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.pdf$
+^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_cache
+^inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_files
 ^vignettes/.build.timestamp$
 ^vignettes/.*_cache$
 ^vignettes/.*cache$
diff --git a/.gitignore b/.gitignore
index 4f479259..9808896d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,8 +1,8 @@
 docs/articles/*_cache/
 install.log
-inst/rmarkdown/templates/hier/skeleton/skeleton_cache
-inst/rmarkdown/templates/hier/skeleton/skeleton_files
-inst/rmarkdown/templates/hier/skeleton/skeleton.pdf
+inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_cache
+inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton_files
+inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.pdf
 mkin*.tar.gz
 mkin.tar
 mkin.Rcheck
diff --git a/GNUmakefile b/GNUmakefile
index 4680514c..4a33c538 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -24,8 +24,8 @@ pkgfiles = \
 	DESCRIPTION \
 	inst/WORDLIST \
 	inst/dataset_generation/* \
-	inst/rmarkdown/templates/hier/template.yaml \
-	inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd \
+	inst/rmarkdown/templates/hierarchical_kinetics/template.yaml \
+	inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd \
 	inst/testdata/* \
 	man/* \
 	NAMESPACE \
diff --git a/NAMESPACE b/NAMESPACE
index 107ffc54..84d6b713 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -96,6 +96,7 @@ export(create_deg_func)
 export(endpoints)
 export(f_time_norm_focus)
 export(get_deg_func)
+export(hierarchical_kinetics)
 export(illparms)
 export(ilr)
 export(intervals)
@@ -187,6 +188,7 @@ importFrom(stats,qf)
 importFrom(stats,qlogis)
 importFrom(stats,qnorm)
 importFrom(stats,qt)
+importFrom(stats,quantile)
 importFrom(stats,residuals)
 importFrom(stats,rnorm)
 importFrom(stats,shapiro.test)
diff --git a/NEWS.md b/NEWS.md
index 7e65204f..2c09df35 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,9 @@
 # mkin 1.2.2
 
+- 'inst/rmarkdown/templates/hier': R markdown template to facilitate the application of hierarchical kinetic models.
+
+- 'inst/testdata/lambda-cyhalothrin_soil_efsa_2014.xlsx': Example spreadsheet for use with 'read_spreadsheet()'.
+
 - 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
 
 - 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted. List correlations of random effects if specified by the user in the covariance model.
diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
new file mode 100644
index 00000000..f7ffb333
--- /dev/null
+++ b/R/hierarchical_kinetics.R
@@ -0,0 +1,39 @@
+#' Hierarchical kinetics template
+#'
+#' R markdown format for setting up hierarchical kinetics based on a template
+#' provided with the mkin package.
+#'
+#' @inheritParams rmarkdown::pdf_document
+#' @param ... Arguments to \code{rmarkdown::pdf_document}
+#'
+#' @return R Markdown output format to pass to
+#'   \code{\link[rmarkdown:render]{render}}
+#'
+#' @examples
+#'
+#' \dontrun{
+#' library(rmarkdown)
+#' draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+#' }
+#'
+#' @export
+hierarchical_kinetics <- function(..., keep_tex = FALSE) {
+
+  if (getRversion() < "4.1.0")
+    stop("You need R with version > 4.1.0 to compile this document")
+
+  if (!requireNamespace("knitr")) stop("Please install the knitr package to use this template")
+  if (!requireNamespace("rmarkdown")) stop("Please install the rmarkdown package to use this template")
+  knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE)
+  knitr::opts_chunk$set(fig.align = "center", fig.pos = "H")
+  options(knitr.kable.NA = "")
+
+  fmt <- rmarkdown::pdf_document(...,
+    keep_tex = keep_tex,
+    toc = TRUE,
+    includes = rmarkdown::includes(in_header = "header.tex"),
+    extra_dependencies = c("float", "listing", "framed")
+  )
+
+  return(fmt)
+}
diff --git a/R/parplot.R b/R/parplot.R
index e9c18947..3da4b51a 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -23,7 +23,7 @@
 #' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
 #' doi: 10.1186/s12859-021-04373-4.
 #' @seealso [multistart]
-#' @importFrom stats median
+#' @importFrom stats median quantile
 #' @export
 parplot <- function(object, ...) {
   UseMethod("parplot")
diff --git a/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd
deleted file mode 100644
index df6fa2f9..00000000
--- a/inst/rmarkdown/templates/hier/skeleton/skeleton.Rmd
+++ /dev/null
@@ -1,186 +0,0 @@
----
-title: "Hierarchical kinetic modelling of degradation data"
-author:
-date: Last change on DD MMM YYYY, last compiled on `r format(Sys.time(),
-  "%e %B %Y")`
-output:
-  pdf_document:
-    extra_dependencies: ["float", "listing"]
-toc: yes
-geometry: margin=2cm
----
-
-
-```{r setup, echo = FALSE, cache = FALSE}
-errmods <- c(const = "constant variance", tc = "two-component error")
-
-knitr::opts_chunk$set(
-  comment = "", tidy = FALSE, cache = TRUE, fig.pos = "H", fig.align = "center"
-)
-options(knitr.kable.NA = "")
-
-# Version requirements
-if (getRversion() < "4.1.0")
-  stop("You need R with version > 4.1.0 to compile this document")
-if ((saemix_version <- packageVersion("saemix")) < "3.1") {
-  warning("Your saemix version is ", saemix_version,
-    ", you should preferably use 3.2 to compile this document")
-}
-if ((mkin_version <- packageVersion("mkin")) < "1.2.2") {
-  stop("Your mkin version is ", mkin_version,
-    ", you need at least 1.2.2 to compile this document")
-}
-```
-
-```{r packages, cache = FALSE, message = FALSE, warning = FALSE, echo = FALSE}
-library(mkin)
-library(saemix)
-library(parallel)
-library(knitr)
-```
-
-```{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")
-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}
-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-pathway, results = "asis", echo = FALSE}
-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}
-parallel::stopCluster(cl)
-sessionInfo()
-```
-
diff --git a/inst/rmarkdown/templates/hier/template.yaml b/inst/rmarkdown/templates/hier/template.yaml
deleted file mode 100644
index d8ab6a4d..00000000
--- a/inst/rmarkdown/templates/hier/template.yaml
+++ /dev/null
@@ -1,3 +0,0 @@
-name: Hierarchical kinetics
-description: Hierarchical kinetic modelling of degradation data
-create_dir: true
diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex
new file mode 100644
index 00000000..a2b7ce83
--- /dev/null
+++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/header.tex
@@ -0,0 +1 @@
+\definecolor{shadecolor}{RGB}{248,248,248}
diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
new file mode 100644
index 00000000..74fae164
--- /dev/null
+++ b/inst/rmarkdown/templates/hierarchical_kinetics/skeleton/skeleton.Rmd
@@ -0,0 +1,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()
+```
+
diff --git a/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml b/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml
new file mode 100644
index 00000000..d8ab6a4d
--- /dev/null
+++ b/inst/rmarkdown/templates/hierarchical_kinetics/template.yaml
@@ -0,0 +1,3 @@
+name: Hierarchical kinetics
+description: Hierarchical kinetic modelling of degradation data
+create_dir: true
diff --git a/log/check.log b/log/check.log
index 42365918..aec61e33 100644
--- a/log/check.log
+++ b/log/check.log
@@ -1,5 +1,5 @@
 * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’
-* using R version 4.2.2 (2022-10-31)
+* using R version 4.2.2 Patched (2022-11-10 r83330)
 * using platform: x86_64-pc-linux-gnu (64-bit)
 * using session charset: UTF-8
 * using options ‘--no-tests --as-cran’
@@ -18,7 +18,7 @@ Maintainer: ‘Johannes Ranke ’
 * checking for portable file names ... OK
 * checking for sufficient/correct file permissions ... OK
 * checking serialization versions ... OK
-* checking whether package ‘mkin’ can be installed ... [11s/11s] OK
+* checking whether package ‘mkin’ can be installed ... OK
 * checking installed package size ... OK
 * checking package directory ... OK
 * checking for future file timestamps ... OK
@@ -41,7 +41,14 @@ Maintainer: ‘Johannes Ranke ’
 * checking S3 generic/method consistency ... OK
 * checking replacement functions ... OK
 * checking foreign function calls ... OK
-* checking R code for possible problems ... [19s/19s] OK
+* checking R code for possible problems ... NOTE
+parplot.multistart.saem.mmkin: no visible global function definition
+  for ‘quantile’
+Undefined global functions or variables:
+  quantile
+Consider adding
+  importFrom("stats", "quantile")
+to your NAMESPACE file.
 * checking Rd files ... OK
 * checking Rd metadata ... OK
 * checking Rd line widths ... OK
@@ -57,7 +64,7 @@ Maintainer: ‘Johannes Ranke ’
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... [24s/24s] OK
+* checking examples ... [11s/11s] OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
@@ -69,5 +76,9 @@ Maintainer: ‘Johannes Ranke ’
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: OK
+Status: 1 NOTE
+See
+  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
+for details.
+
 
diff --git a/man/hierarchical_kinetics.Rd b/man/hierarchical_kinetics.Rd
new file mode 100644
index 00000000..4bb82a4c
--- /dev/null
+++ b/man/hierarchical_kinetics.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/hierarchical_kinetics.R
+\name{hierarchical_kinetics}
+\alias{hierarchical_kinetics}
+\title{Hierarchical kinetics template}
+\usage{
+hierarchical_kinetics(..., keep_tex = FALSE)
+}
+\arguments{
+\item{...}{Arguments to \code{rmarkdown::pdf_document}}
+
+\item{keep_tex}{Keep the intermediate tex file used in the conversion to PDF}
+}
+\value{
+R Markdown output format to pass to
+\code{\link[rmarkdown:render]{render}}
+}
+\description{
+R markdown format for setting up hierarchical kinetics based on a template
+provided with the mkin package.
+}
+\examples{
+
+\dontrun{
+library(rmarkdown)
+draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+}
+
+}
-- 
cgit v1.2.3