<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<meta name="description" content="mkin">
<title>Example evaluations of the dimethenamid data from 2018 • mkin</title>
<script src="../../deps/jquery-3.6.0/jquery-3.6.0.min.js"></script><meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<link href="../../deps/bootstrap-5.2.2/bootstrap.min.css" rel="stylesheet">
<script src="../../deps/bootstrap-5.2.2/bootstrap.bundle.min.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- bootstrap-toc --><script src="https://cdn.jsdelivr.net/gh/afeld/bootstrap-toc@v1.0.1/dist/bootstrap-toc.min.js" integrity="sha256-4veVQbu7//Lk5TSmc7YV48MxtMy98e26cf5MrgZYnwo=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- search --><script src="https://cdnjs.cloudflare.com/ajax/libs/fuse.js/6.4.6/fuse.js" integrity="sha512-zv6Ywkjyktsohkbp9bb45V6tEMoWhzFzXis+LrMehmJZZSys19Yxf1dopHx7WzIKxr5tK2dVcYmaCk2uqdjF4A==" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/autocomplete.js/0.38.0/autocomplete.jquery.min.js" integrity="sha512-GU9ayf+66Xx2TmpxqJpliWbT5PiGYxpaG8rfnBEk1LL8l1KGkRShhngwdXK1UgqhAzWpZHSiYPc09/NwDQIGyg==" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mark.js/8.11.1/mark.min.js" integrity="sha512-5CYOlHXGh6QpOFA/TeTylKLWfB3ftPsde7AnmhuitiTX4K5SqCLBeKro6sPS8ilsz1Q4NRx3v8Ko2IBiszzdww==" crossorigin="anonymous"></script><!-- pkgdown --><script src="../../pkgdown.js"></script><meta property="og:title" content="Example evaluations of the dimethenamid data from 2018">
<meta property="og:description" content="mkin">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body>
<a href="#main" class="visually-hidden-focusable">Skip to contents</a>
<nav class="navbar fixed-top navbar-default navbar-expand-lg bg-light"><div class="container">
<a class="navbar-brand me-2" href="../../index.html">mkin</a>
<small class="nav-text text-default me-auto" data-bs-toggle="tooltip" data-bs-placement="bottom" title="Released version">1.2.6</small>
<button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbar" aria-controls="navbar" aria-expanded="false" aria-label="Toggle navigation">
<span class="navbar-toggler-icon"></span>
</button>
<div id="navbar" class="collapse navbar-collapse ms-3">
<ul class="navbar-nav me-auto">
<li class="nav-item">
<a class="nav-link" href="../../reference/index.html">Reference</a>
</li>
<li class="active nav-item dropdown">
<a href="#" class="nav-link dropdown-toggle" data-bs-toggle="dropdown" role="button" aria-expanded="false" aria-haspopup="true" id="dropdown-articles">Articles</a>
<div class="dropdown-menu" aria-labelledby="dropdown-articles">
<a class="dropdown-item" href="../../articles/mkin.html">Introduction to mkin</a>
<div class="dropdown-divider"></div>
<h6 class="dropdown-header" data-toc-skip>Example evaluations with (generalised) nonlinear least squares</h6>
<a class="dropdown-item" href="../../articles/FOCUS_D.html">Example evaluation of FOCUS Example Dataset D</a>
<a class="dropdown-item" href="../../articles/FOCUS_L.html">Example evaluation of FOCUS Laboratory Data L1 to L3</a>
<a class="dropdown-item" href="../../articles/web_only/FOCUS_Z.html">Example evaluation of FOCUS Example Dataset Z</a>
<div class="dropdown-divider"></div>
<h6 class="dropdown-header" data-toc-skip>Example evaluations with hierarchical models (nonlinear mixed-effects models)</h6>
<a class="dropdown-item" href="../../articles/prebuilt/2022_dmta_parent.html">Testing hierarchical parent degradation kinetics with residue data on dimethenamid and dimethenamid-P</a>
<a class="dropdown-item" href="../../articles/prebuilt/2022_dmta_pathway.html">Testing hierarchical pathway kinetics with residue data on dimethenamid and dimethenamid-P</a>
<a class="dropdown-item" href="../../articles/prebuilt/2023_mesotrione_parent.html">Testing covariate modelling in hierarchical parent degradation kinetics with residue data on mesotrione</a>
<a class="dropdown-item" href="../../articles/prebuilt/2022_cyan_pathway.html">Testing hierarchical pathway kinetics with residue data on cyantraniliprole</a>
<a class="dropdown-item" href="../../articles/web_only/dimethenamid_2018.html">Comparison of saemix and nlme evaluations of dimethenamid data from 2018</a>
<a class="dropdown-item" href="../../articles/web_only/multistart.html">Short demo of the multistart method</a>
<div class="dropdown-divider"></div>
<h6 class="dropdown-header" data-toc-skip>Performance</h6>
<a class="dropdown-item" href="../../articles/web_only/compiled_models.html">Performance benefit by using compiled model definitions in mkin</a>
<a class="dropdown-item" href="../../articles/web_only/benchmarks.html">Benchmark timings for mkin</a>
<a class="dropdown-item" href="../../articles/web_only/saem_benchmarks.html">Benchmark timings for saem.mmkin</a>
<div class="dropdown-divider"></div>
<h6 class="dropdown-header" data-toc-skip>Miscellaneous</h6>
<a class="dropdown-item" href="../../articles/twa.html">Calculation of time weighted average concentrations with mkin</a>
<a class="dropdown-item" href="../../articles/web_only/NAFTA_examples.html">Example evaluation of NAFTA SOP Attachment examples</a>
</div>
</li>
<li class="nav-item">
<a class="nav-link" href="../../news/index.html">News</a>
</li>
</ul>
<form class="form-inline my-2 my-lg-0" role="search">
<input type="search" class="form-control me-sm-2" aria-label="Toggle navigation" name="search-input" data-search-index="../../search.json" id="search-input" placeholder="Search for" autocomplete="off">
</form>
<ul class="navbar-nav">
<li class="nav-item">
<a class="external-link nav-link" href="https://github.com/jranke/mkin/" aria-label="github">
<span class="fab fa fab fa-github fa-lg"></span>
</a>
</li>
</ul>
</div>
</div>
</nav><div class="container template-article">
<div class="row">
<main id="main" class="col-md-9"><div class="page-header">
<img src="" class="logo" alt=""><h1>Example evaluations of the dimethenamid data from 2018</h1>
<h4 data-toc-skip class="author">Johannes
Ranke</h4>
<h4 data-toc-skip class="date">Last change 1 July 2022,
built on 16 Nov 2023</h4>
<small class="dont-index">Source: <a href="https://github.com/jranke/mkin/blob/HEAD/vignettes/web_only/dimethenamid_2018.rmd" class="external-link"><code>vignettes/web_only/dimethenamid_2018.rmd</code></a></small>
<div class="d-none name"><code>dimethenamid_2018.rmd</code></div>
</div>
<p><a href="http://www.jrwb.de" class="external-link">Wissenschaftlicher Berater, Kronacher
Str. 12, 79639 Grenzach-Wyhlen, Germany</a></p>
<div class="section level2">
<h2 id="introduction">Introduction<a class="anchor" aria-label="anchor" href="#introduction"></a>
</h2>
<p>A first analysis of the data analysed here was presented in a recent
journal article on nonlinear mixed-effects models in degradation
kinetics <span class="citation">(Ranke et al. 2021)</span>. That
analysis was based on the <code>nlme</code> package and a development
version of the <code>saemix</code> package that was unpublished at the
time. Meanwhile, version 3.0 of the <code>saemix</code> package is
available from the CRAN repository. Also, it turned out that there was
an error in the handling of the Borstel data in the mkin package at the
time, leading to the duplication of a few data points from that soil.
The dataset in the mkin package has been corrected, and the interface to
<code>saemix</code> in the mkin package has been updated to use the
released version.</p>
<p>This vignette is intended to present an up to date analysis of the
data, using the corrected dataset and released versions of
<code>mkin</code> and <code>saemix</code>.</p>
</div>
<div class="section level2">
<h2 id="data">Data<a class="anchor" aria-label="anchor" href="#data"></a>
</h2>
<p>Residue data forming the basis for the endpoints derived in the
conclusion on the peer review of the pesticide risk assessment of
dimethenamid-P published by the European Food Safety Authority (EFSA) in
2018 <span class="citation">(EFSA 2018)</span> were transcribed from the
risk assessment report <span class="citation">(Rapporteur Member State
Germany, Co-Rapporteur Member State Bulgaria 2018)</span> which can be
downloaded from the Open EFSA repository <a href="https://open.efsa.europa.eu" class="external-link">https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716</a>.</p>
<p>The data are <a href="https://pkgdown.jrwb.de/mkin/reference/dimethenamid_2018.html">available
in the mkin package</a>. The following code (hidden by default, please
use the button to the right to show it) treats the data available for
the racemic mixture dimethenamid (DMTA) and its enantiomer
dimethenamid-P (DMTAP) in the same way, as no difference between their
degradation behaviour was identified in the EU risk assessment. The
observation times of each dataset are multiplied with the corresponding
normalisation factor also available in the dataset, in order to make it
possible to describe all datasets with a single set of parameters.</p>
<p>Also, datasets observed in the same soil are merged, resulting in
dimethenamid (DMTA) data from six soils.</p>
<div class="sourceCode" id="cb1"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://pkgdown.jrwb.de/mkin/">mkin</a></span>, quietly <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span>
<span><span class="va">dmta_ds</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html" class="external-link">lapply</a></span><span class="op">(</span><span class="fl">1</span><span class="op">:</span><span class="fl">7</span>, <span class="kw">function</span><span class="op">(</span><span class="va">i</span><span class="op">)</span> <span class="op">{</span></span>
<span> <span class="va">ds_i</span> <span class="op"><-</span> <span class="va">dimethenamid_2018</span><span class="op">$</span><span class="va">ds</span><span class="op">[[</span><span class="va">i</span><span class="op">]</span><span class="op">]</span><span class="op">$</span><span class="va">data</span></span>
<span> <span class="va">ds_i</span><span class="op">[</span><span class="va">ds_i</span><span class="op">$</span><span class="va">name</span> <span class="op">==</span> <span class="st">"DMTAP"</span>, <span class="st">"name"</span><span class="op">]</span> <span class="op"><-</span> <span class="st">"DMTA"</span></span>
<span> <span class="va">ds_i</span><span class="op">$</span><span class="va">time</span> <span class="op"><-</span> <span class="va">ds_i</span><span class="op">$</span><span class="va">time</span> <span class="op">*</span> <span class="va">dimethenamid_2018</span><span class="op">$</span><span class="va">f_time_norm</span><span class="op">[</span><span class="va">i</span><span class="op">]</span></span>
<span> <span class="va">ds_i</span></span>
<span><span class="op">}</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/names.html" class="external-link">names</a></span><span class="op">(</span><span class="va">dmta_ds</span><span class="op">)</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html" class="external-link">sapply</a></span><span class="op">(</span><span class="va">dimethenamid_2018</span><span class="op">$</span><span class="va">ds</span>, <span class="kw">function</span><span class="op">(</span><span class="va">ds</span><span class="op">)</span> <span class="va">ds</span><span class="op">$</span><span class="va">title</span><span class="op">)</span></span>
<span><span class="va">dmta_ds</span><span class="op">[[</span><span class="st">"Elliot"</span><span class="op">]</span><span class="op">]</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html" class="external-link">rbind</a></span><span class="op">(</span><span class="va">dmta_ds</span><span class="op">[[</span><span class="st">"Elliot 1"</span><span class="op">]</span><span class="op">]</span>, <span class="va">dmta_ds</span><span class="op">[[</span><span class="st">"Elliot 2"</span><span class="op">]</span><span class="op">]</span><span class="op">)</span></span>
<span><span class="va">dmta_ds</span><span class="op">[[</span><span class="st">"Elliot 1"</span><span class="op">]</span><span class="op">]</span> <span class="op"><-</span> <span class="cn">NULL</span></span>
<span><span class="va">dmta_ds</span><span class="op">[[</span><span class="st">"Elliot 2"</span><span class="op">]</span><span class="op">]</span> <span class="op"><-</span> <span class="cn">NULL</span></span></code></pre></div>
</div>
<div class="section level2">
<h2 id="parent-degradation">Parent degradation<a class="anchor" aria-label="anchor" href="#parent-degradation"></a>
</h2>
<p>We evaluate the observed degradation of the parent compound using
simple exponential decline (SFO) and biexponential decline (DFOP), using
constant variance (const) and a two-component variance (tc) as error
models.</p>
<div class="section level3">
<h3 id="separate-evaluations">Separate evaluations<a class="anchor" aria-label="anchor" href="#separate-evaluations"></a>
</h3>
<p>As a first step, to get a visual impression of the fit of the
different models, we do separate evaluations for each soil using the
mmkin function from the mkin package:</p>
<div class="sourceCode" id="cb2"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_mkin_const</span> <span class="op"><-</span> <span class="fu"><a href="../../reference/mmkin.html">mmkin</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"SFO"</span>, <span class="st">"DFOP"</span><span class="op">)</span>, <span class="va">dmta_ds</span>,</span>
<span> error_model <span class="op">=</span> <span class="st">"const"</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span>
<span><span class="va">f_parent_mkin_tc</span> <span class="op"><-</span> <span class="fu"><a href="../../reference/mmkin.html">mmkin</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"SFO"</span>, <span class="st">"DFOP"</span><span class="op">)</span>, <span class="va">dmta_ds</span>,</span>
<span> error_model <span class="op">=</span> <span class="st">"tc"</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p>The plot of the individual SFO fits shown below suggests that at
least in some datasets the degradation slows down towards later time
points, and that the scatter of the residuals error is smaller for
smaller values (panel to the right):</p>
<div class="sourceCode" id="cb3"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="fu"><a href="../../reference/mixed.html">mixed</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_mkin_sfo_const-1.png" width="700"></p>
<p>Using biexponential decline (DFOP) results in a slightly more random
scatter of the residuals:</p>
<div class="sourceCode" id="cb4"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="fu"><a href="../../reference/mixed.html">mixed</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const-1.png" width="700"></p>
<p>The population curve (bold line) in the above plot results from
taking the mean of the individual transformed parameters, i.e. of log k1
and log k2, as well as of the logit of the g parameter of the DFOP
model). Here, this procedure does not result in parameters that
represent the degradation well, because in some datasets the fitted
value for k2 is extremely close to zero, leading to a log k2 value that
dominates the average. This is alleviated if only rate constants that
pass the t-test for significant difference from zero (on the
untransformed scale) are considered in the averaging:</p>
<div class="sourceCode" id="cb5"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="fu"><a href="../../reference/mixed.html">mixed</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span>, test_log_parms <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_const_test-1.png" width="700"></p>
<p>While this is visually much more satisfactory, such an average
procedure could introduce a bias, as not all results from the individual
fits enter the population curve with the same weight. This is where
nonlinear mixed-effects models can help out by treating all datasets
with equally by fitting a parameter distribution model together with the
degradation model and the error model (see below).</p>
<p>The remaining trend of the residuals to be higher for higher
predicted residues is reduced by using the two-component error
model:</p>
<div class="sourceCode" id="cb6"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="fu"><a href="../../reference/mixed.html">mixed</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span>, test_log_parms <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_mkin_dfop_tc_test-1.png" width="700"></p>
<p>However, note that in the case of using this error model, the fits to
the Flaach and BBA 2.3 datasets appear to be ill-defined, indicated by
the fact that they did not converge:</p>
<div class="sourceCode" id="cb7"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span></span></code></pre></div>
<pre><code><mmkin> object
Status of individual fits:
dataset
model Calke Borstel Flaach BBA 2.2 BBA 2.3 Elliot
DFOP OK OK OK OK C OK
C: Optimisation did not converge:
iteration limit reached without convergence (10)
OK: No warnings</code></pre>
</div>
<div class="section level3">
<h3 id="nonlinear-mixed-effects-models">Nonlinear mixed-effects models<a class="anchor" aria-label="anchor" href="#nonlinear-mixed-effects-models"></a>
</h3>
<p>Instead of taking a model selection decision for each of the
individual fits, we fit nonlinear mixed-effects models (using different
fitting algorithms as implemented in different packages) and do model
selection using all available data at the same time. In order to make
sure that these decisions are not unduly influenced by the type of
algorithm used, by implementation details or by the use of wrong control
parameters, we compare the model selection results obtained with
different R packages, with different algorithms and checking control
parameters.</p>
<div class="section level4">
<h4 id="nlme">nlme<a class="anchor" aria-label="anchor" href="#nlme"></a>
</h4>
<p>The nlme package was the first R extension providing facilities to
fit nonlinear mixed-effects models. We would like to do model selection
from all four combinations of degradation models and error models based
on the AIC. However, fitting the DFOP model with constant variance and
using default control parameters results in an error, signalling that
the maximum number of 50 iterations was reached, potentially indicating
overparameterisation. Nevertheless, the algorithm converges when the
two-component error model is used in combination with the DFOP model.
This can be explained by the fact that the smaller residues observed at
later sampling times get more weight when using the two-component error
model which will counteract the tendency of the algorithm to try
parameter combinations unsuitable for fitting these data.</p>
<div class="sourceCode" id="cb9"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://svn.r-project.org/R-packages/trunk/nlme/" class="external-link">nlme</a></span><span class="op">)</span></span>
<span><span class="va">f_parent_nlme_sfo_const</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span><span class="op">)</span></span>
<span><span class="co"># f_parent_nlme_dfop_const <- nlme(f_parent_mkin_const["DFOP", ])</span></span>
<span><span class="va">f_parent_nlme_sfo_tc</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span><span class="op">)</span></span>
<span><span class="va">f_parent_nlme_dfop_tc</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span></span></code></pre></div>
<p>Note that a certain degree of overparameterisation is also indicated
by a warning obtained when fitting DFOP with the two-component error
model (‘false convergence’ in the ‘LME step’ in iteration 3). However,
as this warning does not occur in later iterations, and specifically not
in the last of the 5 iterations, we can ignore this warning.</p>
<p>The model comparison function of the nlme package can directly be
applied to these fits showing a much lower AIC for the DFOP model fitted
with the two-component error model. Also, the likelihood ratio test
indicates that this difference is significant as the p-value is below
0.0001.</p>
<div class="sourceCode" id="cb10"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/stats/anova.html" class="external-link">anova</a></span><span class="op">(</span></span>
<span> <span class="va">f_parent_nlme_sfo_const</span>, <span class="va">f_parent_nlme_sfo_tc</span>, <span class="va">f_parent_nlme_dfop_tc</span></span>
<span><span class="op">)</span></span></code></pre></div>
<pre><code> Model df AIC BIC logLik Test L.Ratio p-value
f_parent_nlme_sfo_const 1 5 796.60 811.82 -393.30
f_parent_nlme_sfo_tc 2 6 798.60 816.86 -393.30 1 vs 2 0.00 0.998
f_parent_nlme_dfop_tc 3 10 671.91 702.34 -325.95 2 vs 3 134.69 <.0001</code></pre>
<p>In addition to these fits, attempts were also made to include
correlations between random effects by using the log Cholesky
parameterisation of the matrix specifying them. The code used for these
attempts can be made visible below.</p>
<div class="sourceCode" id="cb12"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_nlme_sfo_const_logchol</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span>,</span>
<span> random <span class="op">=</span> <span class="fu">nlme</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/nlme/man/pdLogChol.html" class="external-link">pdLogChol</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="va">DMTA_0</span> <span class="op">~</span> <span class="fl">1</span>, <span class="va">log_k_DMTA</span> <span class="op">~</span> <span class="fl">1</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/stats/anova.html" class="external-link">anova</a></span><span class="op">(</span><span class="va">f_parent_nlme_sfo_const</span>, <span class="va">f_parent_nlme_sfo_const_logchol</span><span class="op">)</span></span>
<span><span class="va">f_parent_nlme_sfo_tc_logchol</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span>,</span>
<span> random <span class="op">=</span> <span class="fu">nlme</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/nlme/man/pdLogChol.html" class="external-link">pdLogChol</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="va">DMTA_0</span> <span class="op">~</span> <span class="fl">1</span>, <span class="va">log_k_DMTA</span> <span class="op">~</span> <span class="fl">1</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/stats/anova.html" class="external-link">anova</a></span><span class="op">(</span><span class="va">f_parent_nlme_sfo_tc</span>, <span class="va">f_parent_nlme_sfo_tc_logchol</span><span class="op">)</span></span>
<span><span class="va">f_parent_nlme_dfop_tc_logchol</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/nlme/man/nlme.html" class="external-link">nlme</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span>,</span>
<span> random <span class="op">=</span> <span class="fu">nlme</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/nlme/man/pdLogChol.html" class="external-link">pdLogChol</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="va">DMTA_0</span> <span class="op">~</span> <span class="fl">1</span>, <span class="va">log_k1</span> <span class="op">~</span> <span class="fl">1</span>, <span class="va">log_k2</span> <span class="op">~</span> <span class="fl">1</span>, <span class="va">g_qlogis</span> <span class="op">~</span> <span class="fl">1</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/stats/anova.html" class="external-link">anova</a></span><span class="op">(</span><span class="va">f_parent_nlme_dfop_tc</span>, <span class="va">f_parent_nlme_dfop_tc_logchol</span><span class="op">)</span></span></code></pre></div>
<p>While the SFO variants converge fast, the additional parameters
introduced by this lead to convergence warnings for the DFOP model. The
model comparison clearly show that adding correlations between random
effects does not improve the fits.</p>
<p>The selected model (DFOP with two-component error) fitted to the data
assuming no correlations between random effects is shown below.</p>
<div class="sourceCode" id="cb13"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="va">f_parent_nlme_dfop_tc</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/plot_parent_nlme-1.png" width="700"></p>
</div>
<div class="section level4">
<h4 id="saemix">saemix<a class="anchor" aria-label="anchor" href="#saemix"></a>
</h4>
<p>The saemix package provided the first Open Source implementation of
the Stochastic Approximation to the Expectation Maximisation (SAEM)
algorithm. SAEM fits of degradation models can be conveniently performed
using an interface to the saemix package available in current
development versions of the mkin package.</p>
<p>The corresponding SAEM fits of the four combinations of degradation
and error models are fitted below. As there is no convergence criterion
implemented in the saemix package, the convergence plots need to be
manually checked for every fit. We define control settings that work
well for all the parent data fits shown in this vignette.</p>
<div class="sourceCode" id="cb14"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va">saemix</span><span class="op">)</span></span>
<span><span class="va">saemix_control</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/saemix/man/saemixControl.html" class="external-link">saemixControl</a></span><span class="op">(</span>nbiter.saemix <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">800</span>, <span class="fl">300</span><span class="op">)</span>, nb.chains <span class="op">=</span> <span class="fl">15</span>,</span>
<span> print <span class="op">=</span> <span class="cn">FALSE</span>, save <span class="op">=</span> <span class="cn">FALSE</span>, save.graphs <span class="op">=</span> <span class="cn">FALSE</span>, displayProgress <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span>
<span><span class="va">saemix_control_moreiter</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/saemix/man/saemixControl.html" class="external-link">saemixControl</a></span><span class="op">(</span>nbiter.saemix <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">1600</span>, <span class="fl">300</span><span class="op">)</span>, nb.chains <span class="op">=</span> <span class="fl">15</span>,</span>
<span> print <span class="op">=</span> <span class="cn">FALSE</span>, save <span class="op">=</span> <span class="cn">FALSE</span>, save.graphs <span class="op">=</span> <span class="cn">FALSE</span>, displayProgress <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span>
<span><span class="va">saemix_control_10k</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/pkg/saemix/man/saemixControl.html" class="external-link">saemixControl</a></span><span class="op">(</span>nbiter.saemix <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">10000</span>, <span class="fl">300</span><span class="op">)</span>, nb.chains <span class="op">=</span> <span class="fl">15</span>,</span>
<span> print <span class="op">=</span> <span class="cn">FALSE</span>, save <span class="op">=</span> <span class="cn">FALSE</span>, save.graphs <span class="op">=</span> <span class="cn">FALSE</span>, displayProgress <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span></code></pre></div>
<p>The convergence plot for the SFO model using constant variance is
shown below.</p>
<div class="sourceCode" id="cb15"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_sfo_const</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span>,</span>
<span> control <span class="op">=</span> <span class="va">saemix_control</span>, transformations <span class="op">=</span> <span class="st">"saemix"</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="va">f_parent_saemix_sfo_const</span><span class="op">$</span><span class="va">so</span>, plot.type <span class="op">=</span> <span class="st">"convergence"</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_const-1.png" width="700"></p>
<p>Obviously the selected number of iterations is sufficient to reach
convergence. This can also be said for the SFO fit using the
two-component error model.</p>
<div class="sourceCode" id="cb16"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_sfo_tc</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"SFO"</span>, <span class="op">]</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span>,</span>
<span> control <span class="op">=</span> <span class="va">saemix_control</span>, transformations <span class="op">=</span> <span class="st">"saemix"</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="va">f_parent_saemix_sfo_tc</span><span class="op">$</span><span class="va">so</span>, plot.type <span class="op">=</span> <span class="st">"convergence"</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_saemix_sfo_tc-1.png" width="700"></p>
<p>When fitting the DFOP model with constant variance (see below),
parameter convergence is not as unambiguous.</p>
<div class="sourceCode" id="cb17"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_dfop_const</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_const</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span>,</span>
<span> control <span class="op">=</span> <span class="va">saemix_control</span>, transformations <span class="op">=</span> <span class="st">"saemix"</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_const</span><span class="op">$</span><span class="va">so</span>, plot.type <span class="op">=</span> <span class="st">"convergence"</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_const-1.png" width="700"></p>
<div class="sourceCode" id="cb18"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_const</span><span class="op">)</span></span></code></pre></div>
<pre><code>Kinetic nonlinear mixed-effects model fit by SAEM
Structural model:
d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
* DMTA
Data:
155 observations of 1 variable(s) grouped in 6 datasets
Likelihood computed by importance sampling
AIC BIC logLik
706 704 -344
Fitted parameters:
estimate lower upper
DMTA_0 97.99583 96.50079 99.4909
k1 0.06377 0.03432 0.0932
k2 0.00848 0.00444 0.0125
g 0.95701 0.91313 1.0009
a.1 1.82141 1.60516 2.0377
SD.DMTA_0 1.64787 0.45729 2.8384
SD.k1 0.57439 0.24731 0.9015
SD.k2 0.03296 -2.50524 2.5712
SD.g 1.10266 0.32354 1.8818</code></pre>
<p>While the other parameters converge to credible values, the variance
of k2 (<code>omega2.k2</code>) converges to a very small value. The
printout of the <code>saem.mmkin</code> model shows that the estimated
standard deviation of k2 across the population of soils
(<code>SD.k2</code>) is ill-defined, indicating overparameterisation of
this model.</p>
<p>When the DFOP model is fitted with the two-component error model, we
also observe that the estimated variance of k2 becomes very small, while
being ill-defined, as illustrated by the excessive confidence interval
of <code>SD.k2</code>.</p>
<div class="sourceCode" id="cb20"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_dfop_tc</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span>,</span>
<span> control <span class="op">=</span> <span class="va">saemix_control</span>, transformations <span class="op">=</span> <span class="st">"saemix"</span><span class="op">)</span></span>
<span><span class="va">f_parent_saemix_dfop_tc_moreiter</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span>, quiet <span class="op">=</span> <span class="cn">TRUE</span>,</span>
<span> control <span class="op">=</span> <span class="va">saemix_control_moreiter</span>, transformations <span class="op">=</span> <span class="st">"saemix"</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/plot.html" class="external-link">plot</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span>, plot.type <span class="op">=</span> <span class="st">"convergence"</span><span class="op">)</span></span></code></pre></div>
<p><img src="dimethenamid_2018_files/figure-html/f_parent_saemix_dfop_tc-1.png" width="700"></p>
<div class="sourceCode" id="cb21"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">)</span></span></code></pre></div>
<pre><code>Kinetic nonlinear mixed-effects model fit by SAEM
Structural model:
d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
* DMTA
Data:
155 observations of 1 variable(s) grouped in 6 datasets
Likelihood computed by importance sampling
AIC BIC logLik
666 664 -323
Fitted parameters:
estimate lower upper
DMTA_0 98.24165 96.29190 100.1914
k1 0.06421 0.03352 0.0949
k2 0.00866 0.00617 0.0111
g 0.95340 0.91218 0.9946
a.1 1.06463 0.86503 1.2642
b.1 0.02964 0.02259 0.0367
SD.DMTA_0 2.03611 0.40416 3.6681
SD.k1 0.59534 0.25692 0.9338
SD.k2 0.00042 -73.01372 73.0146
SD.g 1.04234 0.37189 1.7128</code></pre>
<p>Doubling the number of iterations in the first phase of the algorithm
leads to a slightly lower likelihood, and therefore to slightly higher
AIC and BIC values. With even more iterations, the algorithm stops with
an error message. This is related to the variance of k2 approximating
zero and has been submitted as a <a href="https://github.com/saemixdevelopment/saemixextension/issues/29" class="external-link">bug
to the saemix package</a>, as the algorithm does not converge in this
case.</p>
<p>An alternative way to fit DFOP in combination with the two-component
error model is to use the model formulation with transformed parameters
as used per default in mkin. When using this option, convergence is
slower, but eventually the algorithm stops as well with the same error
message.</p>
<p>The four combinations (SFO/const, SFO/tc, DFOP/const and DFOP/tc) and
the version with increased iterations can be compared using the model
comparison function of the saemix package:</p>
<div class="sourceCode" id="cb23"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">AIC_parent_saemix</span> <span class="op"><-</span> <span class="fu">saemix</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/saemix/man/compare.saemix.html" class="external-link">compare.saemix</a></span><span class="op">(</span></span>
<span> <span class="va">f_parent_saemix_sfo_const</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_sfo_tc</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_dfop_const</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_dfop_tc_moreiter</span><span class="op">$</span><span class="va">so</span><span class="op">)</span></span></code></pre></div>
<pre><code>Likelihoods calculated by importance sampling</code></pre>
<div class="sourceCode" id="cb25"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/colnames.html" class="external-link">rownames</a></span><span class="op">(</span><span class="va">AIC_parent_saemix</span><span class="op">)</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span></span>
<span> <span class="st">"SFO const"</span>, <span class="st">"SFO tc"</span>, <span class="st">"DFOP const"</span>, <span class="st">"DFOP tc"</span>, <span class="st">"DFOP tc more iterations"</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">AIC_parent_saemix</span><span class="op">)</span></span></code></pre></div>
<pre><code> AIC BIC
SFO const 796.38 795.34
SFO tc 798.38 797.13
DFOP const 705.75 703.88
DFOP tc 665.67 663.59
DFOP tc more iterations 665.85 663.76</code></pre>
<p>In order to check the influence of the likelihood calculation
algorithms implemented in saemix, the likelihood from Gaussian
quadrature is added to the best fit, and the AIC values obtained from
the three methods are compared.</p>
<div class="sourceCode" id="cb27"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span> <span class="op"><-</span></span>
<span> <span class="fu">saemix</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/saemix/man/llgq.saemix.html" class="external-link">llgq.saemix</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span><span class="op">)</span></span>
<span><span class="va">AIC_parent_saemix_methods</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span></span>
<span> is <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"is"</span><span class="op">)</span>,</span>
<span> gq <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"gq"</span><span class="op">)</span>,</span>
<span> lin <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"lin"</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">AIC_parent_saemix_methods</span><span class="op">)</span></span></code></pre></div>
<pre><code> is gq lin
665.67 665.74 665.13 </code></pre>
<p>The AIC values based on importance sampling and Gaussian quadrature
are very similar. Using linearisation is known to be less accurate, but
still gives a similar value.</p>
<p>In order to illustrate that the comparison of the three method
depends on the degree of convergence obtained in the fit, the same
comparison is shown below for the fit using the defaults for the number
of iterations and the number of MCMC chains.</p>
<p>When using OpenBlas for linear algebra, there is a large difference
in the values obtained with Gaussian quadrature, so the larger number of
iterations makes a lot of difference. When using the LAPACK version
coming with Debian Bullseye, the AIC based on Gaussian quadrature is
almost the same as the one obtained with the other methods, also when
using defaults for the fit.</p>
<div class="sourceCode" id="cb29"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">f_parent_saemix_dfop_tc_defaults</span> <span class="op"><-</span> <span class="fu">mkin</span><span class="fu">::</span><span class="fu"><a href="../../reference/saem.html">saem</a></span><span class="op">(</span><span class="va">f_parent_mkin_tc</span><span class="op">[</span><span class="st">"DFOP"</span>, <span class="op">]</span><span class="op">)</span></span>
<span><span class="va">f_parent_saemix_dfop_tc_defaults</span><span class="op">$</span><span class="va">so</span> <span class="op"><-</span></span>
<span> <span class="fu">saemix</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/saemix/man/llgq.saemix.html" class="external-link">llgq.saemix</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc_defaults</span><span class="op">$</span><span class="va">so</span><span class="op">)</span></span>
<span><span class="va">AIC_parent_saemix_methods_defaults</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span></span>
<span> is <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc_defaults</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"is"</span><span class="op">)</span>,</span>
<span> gq <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc_defaults</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"gq"</span><span class="op">)</span>,</span>
<span> lin <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_saemix_dfop_tc_defaults</span><span class="op">$</span><span class="va">so</span>, method <span class="op">=</span> <span class="st">"lin"</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/print.html" class="external-link">print</a></span><span class="op">(</span><span class="va">AIC_parent_saemix_methods_defaults</span><span class="op">)</span></span></code></pre></div>
<pre><code> is gq lin
670.09 669.37 671.29 </code></pre>
</div>
</div>
<div class="section level3">
<h3 id="comparison">Comparison<a class="anchor" aria-label="anchor" href="#comparison"></a>
</h3>
<p>The following table gives the AIC values obtained with both backend
packages using the same control parameters (800 iterations burn-in, 300
iterations second phase, 15 chains).</p>
<div class="sourceCode" id="cb31"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">AIC_all</span> <span class="op"><-</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html" class="external-link">data.frame</a></span><span class="op">(</span></span>
<span> check.names <span class="op">=</span> <span class="cn">FALSE</span>,</span>
<span> <span class="st">"Degradation model"</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"SFO"</span>, <span class="st">"SFO"</span>, <span class="st">"DFOP"</span>, <span class="st">"DFOP"</span><span class="op">)</span>,</span>
<span> <span class="st">"Error model"</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"const"</span>, <span class="st">"tc"</span>, <span class="st">"const"</span>, <span class="st">"tc"</span><span class="op">)</span>,</span>
<span> nlme <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_nlme_sfo_const</span><span class="op">)</span>, <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_nlme_sfo_tc</span><span class="op">)</span>, <span class="cn">NA</span>, <span class="fu"><a href="https://rdrr.io/r/stats/AIC.html" class="external-link">AIC</a></span><span class="op">(</span><span class="va">f_parent_nlme_dfop_tc</span><span class="op">)</span><span class="op">)</span>,</span>
<span> saemix_lin <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html" class="external-link">sapply</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="va">f_parent_saemix_sfo_const</span><span class="op">$</span><span class="va">so</span>, <span class="va">f_parent_saemix_sfo_tc</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_dfop_const</span><span class="op">$</span><span class="va">so</span>, <span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span><span class="op">)</span>, <span class="va">AIC</span>, method <span class="op">=</span> <span class="st">"lin"</span><span class="op">)</span>,</span>
<span> saemix_is <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html" class="external-link">sapply</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="va">f_parent_saemix_sfo_const</span><span class="op">$</span><span class="va">so</span>, <span class="va">f_parent_saemix_sfo_tc</span><span class="op">$</span><span class="va">so</span>,</span>
<span> <span class="va">f_parent_saemix_dfop_const</span><span class="op">$</span><span class="va">so</span>, <span class="va">f_parent_saemix_dfop_tc</span><span class="op">$</span><span class="va">so</span><span class="op">)</span>, <span class="va">AIC</span>, method <span class="op">=</span> <span class="st">"is"</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/pkg/knitr/man/kable.html" class="external-link">kable</a></span><span class="op">(</span><span class="va">AIC_all</span><span class="op">)</span></span></code></pre></div>
<table class="table">
<thead><tr class="header">
<th align="left">Degradation model</th>
<th align="left">Error model</th>
<th align="right">nlme</th>
<th align="right">saemix_lin</th>
<th align="right">saemix_is</th>
</tr></thead>
<tbody>
<tr class="odd">
<td align="left">SFO</td>
<td align="left">const</td>
<td align="right">796.60</td>
<td align="right">796.60</td>
<td align="right">796.38</td>
</tr>
<tr class="even">
<td align="left">SFO</td>
<td align="left">tc</td>
<td align="right">798.60</td>
<td align="right">798.60</td>
<td align="right">798.38</td>
</tr>
<tr class="odd">
<td align="left">DFOP</td>
<td align="left">const</td>
<td align="right">NA</td>
<td align="right">704.95</td>
<td align="right">705.75</td>
</tr>
<tr class="even">
<td align="left">DFOP</td>
<td align="left">tc</td>
<td align="right">671.91</td>
<td align="right">665.13</td>
<td align="right">665.67</td>
</tr>
</tbody>
</table>
</div>
</div>
<div class="section level2">
<h2 id="conclusion">Conclusion<a class="anchor" aria-label="anchor" href="#conclusion"></a>
</h2>
<p>A more detailed analysis of the dimethenamid dataset confirmed that
the DFOP model provides the most appropriate description of the decline
of the parent compound in these data. On the other hand, closer
inspection of the results revealed that the variability of the k2
parameter across the population of soils is ill-defined. This coincides
with the observation that this parameter cannot robustly be quantified
for some of the soils.</p>
<p>Regarding the regulatory use of these data, it is claimed that an
improved characterisation of the mean parameter values across the
population is obtained using the nonlinear mixed-effects models
presented here. However, attempts to quantify the variability of the
slower rate constant of the biphasic decline of dimethenamid indicate
that the data are not sufficient to characterise this variability to a
satisfactory precision.</p>
</div>
<div class="section level2">
<h2 id="session-info">Session Info<a class="anchor" aria-label="anchor" href="#session-info"></a>
</h2>
<div class="sourceCode" id="cb32"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/sessionInfo.html" class="external-link">sessionInfo</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code>R version 4.3.2 (2023-10-31)
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=C 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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] saemix_3.2 npde_3.3 nlme_3.1-163 mkin_1.2.6 knitr_1.42
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.21-8 digest_0.6.31 magrittr_2.0.3 evaluate_0.21
[9] grid_4.3.2 fastmap_1.1.1 rprojroot_2.0.3 jsonlite_1.8.4
[13] DBI_1.1.3 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.2 parallel_4.3.2
[29] memoise_2.0.1 dplyr_1.1.2 colorspace_2.1-0 ggplot2_3.4.2
[33] vctrs_0.6.2 R6_2.5.1 zoo_1.8-12 lifecycle_1.0.3
[37] stringr_1.5.0 fs_1.6.2 MASS_7.3-60 ragg_1.2.5
[41] pkgconfig_2.0.3 desc_1.4.2 pkgdown_2.0.7 bslib_0.4.2
[45] pillar_1.9.0 gtable_0.3.3 glue_1.6.2 systemfonts_1.0.4
[49] highr_0.10 xfun_0.39 tibble_3.2.1 lmtest_0.9-40
[53] tidyselect_1.2.0 htmltools_0.5.5 rmarkdown_2.21 compiler_4.3.2 </code></pre>
</div>
<div class="section level2">
<h2 id="references">References<a class="anchor" aria-label="anchor" href="#references"></a>
</h2>
<!-- vim: set foldmethod=syntax: -->
<div id="refs" class="references csl-bib-body hanging-indent">
<div id="ref-efsa_2018_dimethenamid" class="csl-entry">
EFSA. 2018. <span>“Peer Review of the Pesticide Risk Assessment of the
Active Substance Dimethenamid-p.”</span> <em>EFSA Journal</em> 16: 5211.
</div>
<div id="ref-ranke2021" class="csl-entry">
Ranke, Johannes, Janina Wöltjen, Jana Schmidt, and Emmanuelle Comets.
2021. <span>“Taking Kinetic Evaluations of Degradation Data to the Next
Level with Nonlinear Mixed-Effects Models.”</span> <em>Environments</em>
8 (8). <a href="https://doi.org/10.3390/environments8080071" class="external-link">https://doi.org/10.3390/environments8080071</a>.
</div>
<div id="ref-dimethenamid_rar_2018_b8" class="csl-entry">
Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria.
2018. <span>“<span class="nocase">Renewal Assessment Report
Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour, Rev. 2 -
November 2017</span>.”</span> <a href="https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716" class="external-link">https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716</a>.
</div>
</div>
</div>
</main><aside class="col-md-3"><nav id="toc"><h2>On this page</h2>
</nav></aside>
</div>
<footer><div class="pkgdown-footer-left">
<p></p>
<p>Developed by Johannes Ranke.</p>
</div>
<div class="pkgdown-footer-right">
<p></p>
<p>Site built with <a href="https://pkgdown.r-lib.org/" class="external-link">pkgdown</a> 2.0.7.</p>
</div>
</footer>
</div>
</body>
</html>
