From 047d048b89e167fb354b45cd7c6b719b9f4cdd28 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 17 Sep 2021 08:47:09 +0200 Subject: Put the AIC comparison in a subsubsection --- docs/dev/articles/web_only/dimethenamid_2018.html | 46 ++++++++++++----------- docs/dev/pkgdown.yml | 2 +- vignettes/web_only/dimethenamid_2018.html | 23 +++++++----- vignettes/web_only/dimethenamid_2018.rmd | 20 +++++++--- 4 files changed, 54 insertions(+), 37 deletions(-) diff --git a/docs/dev/articles/web_only/dimethenamid_2018.html b/docs/dev/articles/web_only/dimethenamid_2018.html index b35d8210..26b352e1 100644 --- a/docs/dev/articles/web_only/dimethenamid_2018.html +++ b/docs/dev/articles/web_only/dimethenamid_2018.html @@ -101,7 +101,7 @@

Example evaluations of the dimethenamid data from 2018

Johannes Ranke

-

Last change 16 September 2021, built on 16 Sep 2021

+

Last change 17 September 2021, built on 17 Sep 2021

Source: vignettes/web_only/dimethenamid_2018.rmd @@ -114,7 +114,7 @@

Introduction

-

During the preparation of the journal article on nonlinear mixed-effects models in degradation kinetics (Ranke et al. 2021) and the analysis of the dimethenamid degradation data analysed therein, a need for a more detailed analysis using not only nlme and saemix, but also nlmixr for fitting the mixed-effects models was identified.

+

During the preparation of the journal article on nonlinear mixed-effects models in degradation kinetics (Ranke et al. 2021) and the analysis of the dimethenamid degradation data analysed therein, a need for a more detailed analysis using not only nlme and saemix, but also nlmixr for fitting the mixed-effects models was identified, as many model variants do not converge when fitted with nlme, and not all relevant error models can be fitted with saemix.

This vignette is an attempt to satisfy this need.

@@ -124,7 +124,7 @@

The data are available in the mkin package. 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.

Also, datasets observed in the same soil are merged, resulting in dimethenamid (DMTA) data from six soils.

-library(mkin)
+library(mkin, quietly = TRUE)
 dmta_ds <- lapply(1:7, function(i) {
   ds_i <- dimethenamid_2018$ds[[i]]$data
   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
@@ -258,14 +258,14 @@ f_parent_nlme_dfop_tc       3 10 671.91 702.34 -325.96 2 vs 3  134.69  <.0001
 
 f_parent_saemix_dfop_tc$so <-
   llgq.saemix(f_parent_saemix_dfop_tc$so)
-AIC(f_parent_saemix_dfop_tc$so)

-
[1] 666.1

-
-AIC(f_parent_saemix_dfop_tc$so, method = "gq")

-
[1] 666.03

-
-AIC(f_parent_saemix_dfop_tc$so, method = "lin")

-
[1] 665.48

+AIC_parent_saemix_methods <- c(
+  is = AIC(f_parent_saemix_dfop_tc$so, method = "is"),
+  gq = AIC(f_parent_saemix_dfop_tc$so, method = "gq"),
+  lin = AIC(f_parent_saemix_dfop_tc$so, method = "lin")
+)
+print(AIC_parent_saemix_methods)
+
    is     gq    lin 
+666.10 666.03 665.48

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.

@@ -273,19 +273,19 @@ f_parent_nlme_dfop_tc 3 10 671.91 702.34 -325.96 2 vs 3 134.69 <.0001 nlmixr

In the last years, a lot of effort has been put into the nlmixr package which is designed for pharmacokinetics, where nonlinear mixed-effects models are routinely used, but which can also be used for related data like chemical degradation data. A current development branch of the mkin package provides an interface between mkin and nlmixr. Here, we check if we get equivalent results when using a refined version of the First Order Conditional Estimation (FOCE) algorithm used in nlme, namely the First Order Conditional Estimation with Interaction (FOCEI), and the SAEM algorithm as implemented in nlmixr.

First, the focei algorithm is used for the four model combinations. A number of warnings are produced with unclear significance.

-
+
 library(nlmixr)
 f_parent_nlmixr_focei_sfo_const <- nlmixr(f_parent_mkin_const["SFO", ], est = "focei")
 f_parent_nlmixr_focei_sfo_tc <- nlmixr(f_parent_mkin_tc["SFO", ], est = "focei")
 f_parent_nlmixr_focei_dfop_const <- nlmixr(f_parent_mkin_const["DFOP", ], est = "focei")
 f_parent_nlmixr_focei_dfop_tc<- nlmixr(f_parent_mkin_tc["DFOP", ], est = "focei")

-
+
 aic_nlmixr_focei <- sapply(
   list(f_parent_nlmixr_focei_sfo_const$nm, f_parent_nlmixr_focei_sfo_tc$nm,
     f_parent_nlmixr_focei_dfop_const$nm, f_parent_nlmixr_focei_dfop_tc$nm),
   AIC)

 
The AIC values are very close to the ones obtained with nlme which are repeated below for convenience.

-
+
 aic_nlme <- sapply(
   list(f_parent_nlme_sfo_const, NA, f_parent_nlme_sfo_tc, f_parent_nlme_dfop_tc),
   function(x) if (is.na(x[1])) NA else AIC(x))
@@ -297,35 +297,35 @@ f_parent_nlme_dfop_tc       3 10 671.91 702.34 -325.96 2 vs 3  134.69  <.0001
   check.names = FALSE
 )

 
Secondly, we use the SAEM estimation routine and check the convergence plots. The control parameters also used for the saemix fits are defined beforehand.

-
+
 nlmixr_saem_control <- saemControl(logLik = TRUE,
   nBurn = 1000, nEm = 300, nmc = 15)

 
The we fit SFO with constant variance

-
+
 f_parent_nlmixr_saem_sfo_const <- nlmixr(f_parent_mkin_const["SFO", ], est = "saem",
   control = nlmixr_saem_control)
 traceplot(f_parent_nlmixr_saem_sfo_const$nm)

 

 
and SFO with two-component error.

-
+
 f_parent_nlmixr_saem_sfo_tc <- nlmixr(f_parent_mkin_tc["SFO", ], est = "saem",
   control = nlmixr_saem_control)
 traceplot(f_parent_nlmixr_saem_sfo_tc$nm)

 

 
For DFOP with constant variance, the convergence plots show considerable instability of the fit, which indicates overparameterisation which was already observed earlier for this model combination.

-
+
 f_parent_nlmixr_saem_dfop_const <- nlmixr(f_parent_mkin_const["DFOP", ], est = "saem",
   control = nlmixr_saem_control)
 traceplot(f_parent_nlmixr_saem_dfop_const$nm)

 

 
For DFOP with two-component error, a less erratic convergence is seen.

-
+
 f_parent_nlmixr_saem_dfop_tc <- nlmixr(f_parent_mkin_tc["DFOP", ], est = "saem",
   control = nlmixr_saem_control)
 traceplot(f_parent_nlmixr_saem_dfop_tc$nm)

 

 
The AIC values are internally calculated using Gaussian quadrature. For an unknown reason, the AIC value obtained for the DFOP fit using constant error is given as Infinity.

-
+
 AIC(f_parent_nlmixr_saem_sfo_const$nm, f_parent_nlmixr_saem_sfo_tc$nm,
   f_parent_nlmixr_saem_dfop_const$nm, f_parent_nlmixr_saem_dfop_tc$nm)

 
                                   df    AIC
@@ -333,8 +333,12 @@ f_parent_nlmixr_saem_sfo_const$nm   5 798.68
 f_parent_nlmixr_saem_sfo_tc$nm      6 808.88
 f_parent_nlmixr_saem_dfop_const$nm  9 815.95
 f_parent_nlmixr_saem_dfop_tc$nm    10 669.57

+

+

+

+Comparison

 
The following table gives the AIC values obtained with the three packages.

-
+
 AIC_all <- data.frame(
   check.names = FALSE,
   "Degradation model" = c("SFO", "SFO", "DFOP", "DFOP"),
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index e932afd0..256387c1 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -11,7 +11,7 @@ articles:
   web_only/benchmarks: benchmarks.html
   web_only/compiled_models: compiled_models.html
   web_only/dimethenamid_2018: dimethenamid_2018.html
-last_built: 2021-09-16T15:10Z
+last_built: 2021-09-17T06:44Z
 urls:
   reference: https://pkgdown.jrwb.de/mkin/reference
   article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/vignettes/web_only/dimethenamid_2018.html b/vignettes/web_only/dimethenamid_2018.html
index aa1905da..a92a720d 100644
--- a/vignettes/web_only/dimethenamid_2018.html
+++ b/vignettes/web_only/dimethenamid_2018.html
@@ -1591,7 +1591,7 @@ div.tocify {
 
 
Example evaluations of the dimethenamid data from 2018

 
Johannes Ranke

-
Last change 16 September 2021, built on 16 Sep 2021

+
Last change 17 September 2021, built on 17 Sep 2021

 
 

 
@@ -1599,7 +1599,7 @@ div.tocify {
 
Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany

 

 
Introduction

-
During the preparation of the journal article on nonlinear mixed-effects models in degradation kinetics (Ranke et al. 2021) and the analysis of the dimethenamid degradation data analysed therein, a need for a more detailed analysis using not only nlme and saemix, but also nlmixr for fitting the mixed-effects models was identified.

+
During the preparation of the journal article on nonlinear mixed-effects models in degradation kinetics (Ranke et al. 2021) and the analysis of the dimethenamid degradation data analysed therein, a need for a more detailed analysis using not only nlme and saemix, but also nlmixr for fitting the mixed-effects models was identified, as many model variants do not converge when fitted with nlme, and not all relevant error models can be fitted with saemix.

 
This vignette is an attempt to satisfy this need.

 

 

@@ -1607,7 +1607,7 @@ div.tocify {
 
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 (EFSA 2018) were transcribed from the risk assessment report (Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria 2018) which can be downloaded from the Open EFSA repository https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716.

 
The data are available in the mkin package. 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.

 
Also, datasets observed in the same soil are merged, resulting in dimethenamid (DMTA) data from six soils.

-
library(mkin)
+
library(mkin, quietly = TRUE)
 dmta_ds <- lapply(1:7, function(i) {
   ds_i <- dimethenamid_2018$ds[[i]]$data
   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
@@ -1720,12 +1720,14 @@ plot(f_parent_saemix_dfop_tc$so, plot.type = "convergence")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.

 
f_parent_saemix_dfop_tc$so <-
   llgq.saemix(f_parent_saemix_dfop_tc$so)
-AIC(f_parent_saemix_dfop_tc$so)

-
[1] 666.1

-
AIC(f_parent_saemix_dfop_tc$so, method = "gq")

-
[1] 666.03

-
AIC(f_parent_saemix_dfop_tc$so, method = "lin")

-
[1] 665.48

+AIC_parent_saemix_methods <- c(
+  is = AIC(f_parent_saemix_dfop_tc$so, method = "is"),
+  gq = AIC(f_parent_saemix_dfop_tc$so, method = "gq"),
+  lin = AIC(f_parent_saemix_dfop_tc$so, method = "lin")
+)
+print(AIC_parent_saemix_methods)

+
    is     gq    lin 
+666.10 666.03 665.48 

 
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.

 

 

@@ -1783,6 +1785,9 @@ f_parent_nlmixr_saem_sfo_const$nm   5 798.68
 f_parent_nlmixr_saem_sfo_tc$nm      6 808.88
 f_parent_nlmixr_saem_dfop_const$nm  9 815.95
 f_parent_nlmixr_saem_dfop_tc$nm    10 669.57

+

+

+
Comparison

 
The following table gives the AIC values obtained with the three packages.

 
AIC_all <- data.frame(
   check.names = FALSE,
diff --git a/vignettes/web_only/dimethenamid_2018.rmd b/vignettes/web_only/dimethenamid_2018.rmd
index e5c8764d..7679edc4 100644
--- a/vignettes/web_only/dimethenamid_2018.rmd
+++ b/vignettes/web_only/dimethenamid_2018.rmd
@@ -1,7 +1,7 @@
 ---
 title: Example evaluations of the dimethenamid data from 2018
 author: Johannes Ranke
-date: Last change 16 September 2021, built on `r format(Sys.Date(), format = "%d %b %Y")`
+date: Last change 17 September 2021, built on `r format(Sys.Date(), format = "%d %b %Y")`
 output:
   html_document:
     toc: true
@@ -32,7 +32,8 @@ During the preparation of the journal article on nonlinear mixed-effects models
 in degradation kinetics [@ranke2021] and the analysis of the dimethenamid
 degradation data analysed therein, a need for a more detailed analysis using
 not only nlme and saemix, but also nlmixr for fitting the mixed-effects models
-was identified.
+was identified, as many model variants do not converge when fitted with nlme,
+and not all relevant error models can be fitted with saemix.
 
 This vignette is an attempt to satisfy this need.
 
@@ -59,7 +60,7 @@ Also, datasets observed in the same soil are merged, resulting in dimethenamid
 (DMTA) data from six soils.
 
 ```{r dimethenamid_data}
-library(mkin)
+library(mkin, quietly = TRUE)
 dmta_ds <- lapply(1:7, function(i) {
   ds_i <- dimethenamid_2018$ds[[i]]$data
   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
@@ -294,9 +295,12 @@ are compared.
 ```{r AIC_parent_saemix_methods, cache = FALSE}
 f_parent_saemix_dfop_tc$so <-
   llgq.saemix(f_parent_saemix_dfop_tc$so)
-AIC(f_parent_saemix_dfop_tc$so)
-AIC(f_parent_saemix_dfop_tc$so, method = "gq")
-AIC(f_parent_saemix_dfop_tc$so, method = "lin")
+AIC_parent_saemix_methods <- c(
+  is = AIC(f_parent_saemix_dfop_tc$so, method = "is"),
+  gq = AIC(f_parent_saemix_dfop_tc$so, method = "gq"),
+  lin = AIC(f_parent_saemix_dfop_tc$so, method = "lin")
+)
+print(AIC_parent_saemix_methods)
 ```
 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
@@ -398,6 +402,8 @@ AIC(f_parent_nlmixr_saem_sfo_const$nm, f_parent_nlmixr_saem_sfo_tc$nm,
   f_parent_nlmixr_saem_dfop_const$nm, f_parent_nlmixr_saem_dfop_tc$nm)
 ```
 
+### Comparison
+
 The following table gives the AIC values obtained with the three packages.
 
 ```{r AIC_all, cache = FALSE}
@@ -416,6 +422,8 @@ AIC_all <- data.frame(
 kable(AIC_all)
 ```
 
+
+
 # References
 
 
-- 
cgit v1.2.1