From e48c1f2ef990a622722e416c8d301430db4f5081 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 16 Mar 2022 14:42:30 +0100 Subject: Use saemix in help page, typo in vignette --- NEWS.md | 6 +- R/dimethenamid_2018.R | 71 ++++----- README.html | 3 +- README.md | 4 + docs/articles/web_only/dimethenamid_2018.html | 4 +- docs/index.html | 1 + docs/news/index.html | 8 +- docs/pkgdown.yml | 2 +- docs/reference/Rplot001.png | Bin 1011 -> 90922 bytes docs/reference/Rplot002.png | Bin 57977 -> 90653 bytes docs/reference/Rplot003.png | Bin 59128 -> 20189 bytes docs/reference/Rplot004.png | Bin 59260 -> 19212 bytes docs/reference/dimethenamid_2018-1.png | Bin 264424 -> 254862 bytes docs/reference/dimethenamid_2018-2.png | Bin 243268 -> 254236 bytes docs/reference/dimethenamid_2018-3.png | Bin 0 -> 7629 bytes docs/reference/dimethenamid_2018.html | 203 ++++++++++++++++++-------- man/dimethenamid_2018.Rd | 71 ++++----- vignettes/web_only/dimethenamid_2018.rmd | 2 +- 18 files changed, 239 insertions(+), 136 deletions(-) create mode 100644 docs/reference/dimethenamid_2018-3.png diff --git a/NEWS.md b/NEWS.md index c0bade8f..e675fc72 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,12 @@ +# mkin 1.1.1 + +- 'dimethenamid_2018': Update example code to use saemix + # mkin 1.1.0 (2022-03-14) ## Mixed-effects models -- Reintroduce the interface to saemix (now on CRAN), in particular the generic function 'saem' with a generator 'saem.mmkin', currently using 'saemix_model' and 'saemix_data', summary and plot methods +- Reintroduce the interface to saemix version 3.0 (now on CRAN), in particular the generic function 'saem' with a generator 'saem.mmkin', currently using 'saemix_model' and 'saemix_data', summary and plot methods - 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for the different mixed-effects model backends diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R index 2bf5fb77..2fdd1981 100644 --- a/R/dimethenamid_2018.R +++ b/R/dimethenamid_2018.R @@ -29,41 +29,42 @@ #' dmta_ds[["Elliot 1"]] <- NULL #' dmta_ds[["Elliot 2"]] <- NULL #' \dontrun{ -#' dfop_sfo3_plus <- mkinmod( -#' DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), -#' M23 = mkinsub("SFO"), -#' M27 = mkinsub("SFO"), -#' M31 = mkinsub("SFO", "M27", sink = FALSE), -#' quiet = TRUE +#' # We don't use DFOP for the parent compound, as this gives numerical +#' # instabilities in the fits +#' sfo_sfo3p <- mkinmod( +#' DMTA = mkinsub("SFO", c("M23", "M27", "M31")), +#' M23 = mkinsub("SFO"), +#' M27 = mkinsub("SFO"), +#' M31 = mkinsub("SFO", "M27", sink = FALSE), +#' quiet = TRUE #' ) -#' f_dmta_mkin_tc <- mmkin( -#' list("DFOP-SFO3+" = dfop_sfo3_plus), -#' dmta_ds, quiet = TRUE, error_model = "tc") -#' nlmixr_model(f_dmta_mkin_tc) -#' # The focei fit takes about four minutes on my system -#' system.time( -#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", -#' control = nlmixr::foceiControl(print = 500)) -#' ) -#' summary(f_dmta_nlmixr_focei) -#' plot(f_dmta_nlmixr_focei) -#' # Using saemix takes about 18 minutes -#' system.time( -#' f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) -#' ) -#' -#' # nlmixr with est = "saem" is pretty fast with default iteration numbers, most -#' # of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end -#' # The likelihood calculated for the nlmixr fit is much lower than that found by saemix -#' # Also, the trace plot and the plot of the individual predictions is not -#' # convincing for the parent. It seems we are fitting an overparameterised -#' # model, so the result we get strongly depends on starting parameters and control settings. -#' system.time( -#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", -#' control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)) -#' ) -#' traceplot(f_dmta_nlmixr_saem$nm) -#' summary(f_dmta_nlmixr_saem) -#' plot(f_dmta_nlmixr_saem) +#' dmta_sfo_sfo3p_tc <- mmkin(list("SFO-SFO3+" = sfo_sfo3p), +#' dmta_ds, error_model = "tc", quiet = TRUE) +#' print(dmta_sfo_sfo3p_tc) +#' # The default (test_log_parms = FALSE) gives an undue +#' # influence of ill-defined rate constants that have +#' # extremely small values: +#' plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = FALSE) +#' # If we disregards ill-defined rate constants, the results +#' # look more plausible, but the truth is likely to be in +#' # between these variants +#' plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = TRUE) +#' # Therefore we use nonlinear mixed-effects models +#' # f_dmta_nlme_tc <- nlme(dmta_sfo_sfo3p_tc) +#' # nlme reaches maxIter = 50 without convergence +#' f_dmta_saem_tc <- saem(dmta_sfo_sfo3p_tc) +#' # I am commenting out the convergence plot as rendering them +#' # with pkgdown fails (at least without further tweaks to the +#' # graphics device used) +#' #saemix::plot(f_dmta_saem_tc$so, plot.type = "convergence") +#' summary(f_dmta_saem_tc) +#' # As the confidence interval for the random effects of DMTA_0 +#' # includes zero, we could try an alternative model without +#' # such random effects +#' # f_dmta_saem_tc_2 <- saem(dmta_sfo_sfo3p_tc, +#' # covariance.model = diag(c(0, rep(1, 7)))) +#' # saemix::plot(f_dmta_saem_tc_2$so, plot.type = "convergence") +#' # This does not perform better judged by AIC and BIC +#' saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so) #' } "dimethenamid_2018" diff --git a/README.html b/README.html index dc5bbcb9..cdd68210 100644 --- a/README.html +++ b/README.html @@ -367,7 +367,7 @@ pre code {

mkin

-

Build Status codecov

+

Build Status codecov

The R package mkin provides calculation routines for the analysis of chemical degradation data, including multicompartment kinetics as needed for modelling the formation and decline of transformation products, or if several degradation compartments are involved.

Installation

@@ -448,6 +448,7 @@ pre code {
  • Project Number 120667 (Development of objective criteria for the evaluation of the visual fit in the kinetic evaluation of degradation data, 2019-2020)
  • Project Number 146839 (Checking the feasibility of using mixed-effects models for the derivation of kinetic modelling parameters from degradation studies, 2020-2021)
    • +

    Thanks are due also to Emmanuelle Comets, maintainer of the saemix package, for the nice collaboration on using the SAEM algorithm and its implementation in saemix for the evaluation of chemical degradation data.

    References

    diff --git a/README.md b/README.md index b864fe79..7c29696a 100644 --- a/README.md +++ b/README.md @@ -212,6 +212,10 @@ to ModelMaker 4.0, 2014-2015) - Project Number 146839 (Checking the feasibility of using mixed-effects models for the derivation of kinetic modelling parameters from degradation studies, 2020-2021) +Thanks are due also to Emmanuelle Comets, maintainer of the saemix package, for +the nice collaboration on using the SAEM algorithm and its implementation in +saemix for the evaluation of chemical degradation data. + ## References diff --git a/docs/articles/web_only/dimethenamid_2018.html b/docs/articles/web_only/dimethenamid_2018.html index 0ba9d5a8..09aa5150 100644 --- a/docs/articles/web_only/dimethenamid_2018.html +++ b/docs/articles/web_only/dimethenamid_2018.html @@ -105,7 +105,7 @@

    Example evaluations of the dimethenamid data from 2018

    Johannes Ranke

    -

    Last change 7 March 2022, built on 07 Mar 2022

    +

    Last change 7 March 2022, built on 16 Mar 2022

    Source: vignettes/web_only/dimethenamid_2018.rmd @@ -426,7 +426,7 @@ DFOP tc more iterations 665.88 663.80

    Conclusion

    -

    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 in some for some of the soils.

    +

    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.

    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.

    diff --git a/docs/index.html b/docs/index.html index 9b5b34a1..4491b50f 100644 --- a/docs/index.html +++ b/docs/index.html @@ -208,6 +208,7 @@
  • Project Number 120667 (Development of objective criteria for the evaluation of the visual fit in the kinetic evaluation of degradation data, 2019-2020)
  • Project Number 146839 (Checking the feasibility of using mixed-effects models for the derivation of kinetic modelling parameters from degradation studies, 2020-2021)
    • +

    Thanks are due also to Emmanuelle Comets, maintainer of the saemix package, for the nice collaboration on using the SAEM algorithm and its implementation in saemix for the evaluation of chemical degradation data.

    References diff --git a/docs/news/index.html b/docs/news/index.html index 075b4e22..13b82597 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -82,10 +82,14 @@

    -

    mkin 1.1.0 (unreleased)

    +

    mkin 1.1.1

    +
    • ‘dimethenamid_2018’: Update example code to use saemix
      • +
    +
    +

    mkin 1.1.0 (2022-03-14)

    Mixed-effects models

    -

    • Reintroduce the interface to saemix (now on CRAN), in particular the generic function ‘saem’ with a generator ‘saem.mmkin’, currently using ‘saemix_model’ and ‘saemix_data’, summary and plot methods

      • +

      • Reintroduce the interface to saemix version 3.0 (now on CRAN), in particular the generic function ‘saem’ with a generator ‘saem.mmkin’, currently using ‘saemix_model’ and ‘saemix_data’, summary and plot methods

      • ‘mean_degparms’: New argument ‘test_log_parms’ that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for the different mixed-effects model backends

      • ‘plot.mixed.mmkin’: Gains arguments ‘test_log_parms’ and ‘conf.level’

      • ‘vignettes/web_only/dimethenamid_2018.rmd’: Example evaluations of the dimethenamid data.

        • diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 59854491..3e26de5b 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -11,7 +11,7 @@ articles: benchmarks: web_only/benchmarks.html compiled_models: web_only/compiled_models.html dimethenamid_2018: web_only/dimethenamid_2018.html -last_built: 2022-03-14T07:56Z +last_built: 2022-03-16T13:42Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png index 17a35806..26efcdd2 100644 Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index ddd3a0d7..065f21f6 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/Rplot003.png b/docs/reference/Rplot003.png index fa29fc43..50a52a1c 100644 Binary files a/docs/reference/Rplot003.png and b/docs/reference/Rplot003.png differ diff --git a/docs/reference/Rplot004.png b/docs/reference/Rplot004.png index 69fdb09d..0f3d8214 100644 Binary files a/docs/reference/Rplot004.png and b/docs/reference/Rplot004.png differ diff --git a/docs/reference/dimethenamid_2018-1.png b/docs/reference/dimethenamid_2018-1.png index 338076a5..be59fc90 100644 Binary files a/docs/reference/dimethenamid_2018-1.png and b/docs/reference/dimethenamid_2018-1.png differ diff --git a/docs/reference/dimethenamid_2018-2.png b/docs/reference/dimethenamid_2018-2.png index 1f41a3a8..0bbb14d8 100644 Binary files a/docs/reference/dimethenamid_2018-2.png and b/docs/reference/dimethenamid_2018-2.png differ diff --git a/docs/reference/dimethenamid_2018-3.png b/docs/reference/dimethenamid_2018-3.png new file mode 100644 index 00000000..71087e85 Binary files /dev/null and b/docs/reference/dimethenamid_2018-3.png differ diff --git a/docs/reference/dimethenamid_2018.html b/docs/reference/dimethenamid_2018.html index f0bc23ee..85e954ac 100644 --- a/docs/reference/dimethenamid_2018.html +++ b/docs/reference/dimethenamid_2018.html @@ -156,65 +156,152 @@ specific pieces of information in the comments.

        dmta_ds[["Elliot 1"]] <- NULL dmta_ds[["Elliot 2"]] <- NULL # \dontrun{ -dfop_sfo3_plus <- mkinmod( - DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), - M23 = mkinsub("SFO"), - M27 = mkinsub("SFO"), - M31 = mkinsub("SFO", "M27", sink = FALSE), - quiet = TRUE +# We don't use DFOP for the parent compound, as this gives numerical +# instabilities in the fits +sfo_sfo3p <- mkinmod( + DMTA = mkinsub("SFO", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + quiet = TRUE ) -f_dmta_mkin_tc <- mmkin( - list("DFOP-SFO3+" = dfop_sfo3_plus), - dmta_ds, quiet = TRUE, error_model = "tc") -nlmixr_model(f_dmta_mkin_tc) -#> Error in nlmixr_model(f_dmta_mkin_tc): could not find function "nlmixr_model" -# The focei fit takes about four minutes on my system -system.time( - f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", - control = nlmixr::foceiControl(print = 500)) -) -#> Error in nlmixr(f_dmta_mkin_tc, est = "focei", control = nlmixr::foceiControl(print = 500)): could not find function "nlmixr" -#> Timing stopped at: 0 0 0 -summary(f_dmta_nlmixr_focei) -#> Error in summary(f_dmta_nlmixr_focei): object 'f_dmta_nlmixr_focei' not found -plot(f_dmta_nlmixr_focei) -#> Error in plot(f_dmta_nlmixr_focei): object 'f_dmta_nlmixr_focei' not found -# Using saemix takes about 18 minutes -system.time( - f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) -) -#> DINTDY- T (=R1) illegal -#> In above message, R1 = 115.507 -#> -#> T not in interval TCUR - HU (= R1) to TCUR (=R2) -#> In above message, R1 = 112.133, R2 = 113.577 -#> -#> DLSODA- At T (=R1), too much accuracy requested -#> for precision of machine.. See TOLSF (=R2) -#> In above message, R1 = 55.3899, R2 = nan -#> -#> Error in out[available, var]: (subscript) logical subscript too long -#> Timing stopped at: 12.76 3.069 11.79 -#> Timing stopped at: 13.77 4.719 12.37 - -# nlmixr with est = "saem" is pretty fast with default iteration numbers, most -# of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end -# The likelihood calculated for the nlmixr fit is much lower than that found by saemix -# Also, the trace plot and the plot of the individual predictions is not -# convincing for the parent. It seems we are fitting an overparameterised -# model, so the result we get strongly depends on starting parameters and control settings. -system.time( - f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", - control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)) -) -#> Error in nlmixr(f_dmta_mkin_tc, est = "saem", control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)): could not find function "nlmixr" -#> Timing stopped at: 0 0 0.001 -traceplot(f_dmta_nlmixr_saem$nm) -#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot" -summary(f_dmta_nlmixr_saem) -#> Error in summary(f_dmta_nlmixr_saem): object 'f_dmta_nlmixr_saem' not found -plot(f_dmta_nlmixr_saem) -#> Error in plot(f_dmta_nlmixr_saem): object 'f_dmta_nlmixr_saem' not found +dmta_sfo_sfo3p_tc <- mmkin(list("SFO-SFO3+" = sfo_sfo3p), + dmta_ds, error_model = "tc", quiet = TRUE) +print(dmta_sfo_sfo3p_tc) +#> <mmkin> object +#> Status of individual fits: +#> +#> dataset +#> model Calke Borstel Flaach BBA 2.2 BBA 2.3 Elliot +#> SFO-SFO3+ OK OK OK OK OK OK +#> +#> OK: No warnings +# The default (test_log_parms = FALSE) gives an undue +# influence of ill-defined rate constants that have +# extremely small values: +plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = FALSE) + +# If we disregards ill-defined rate constants, the results +# look more plausible, but the truth is likely to be in +# between these variants +plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = TRUE) + +# Therefore we use nonlinear mixed-effects models +# f_dmta_nlme_tc <- nlme(dmta_sfo_sfo3p_tc) +# nlme reaches maxIter = 50 without convergence +f_dmta_saem_tc <- saem(dmta_sfo_sfo3p_tc) +# I am commenting out the convergence plot as rendering them +# with pkgdown fails (at least without further tweaks to the +# graphics device used) +#saemix::plot(f_dmta_saem_tc$so, plot.type = "convergence") +summary(f_dmta_saem_tc) +#> saemix version used for fitting: 3.0 +#> mkin version used for pre-fitting: 1.1.0 +#> R version used for fitting: 4.1.3 +#> Date of fit: Wed Mar 16 14:32:04 2022 +#> Date of summary: Wed Mar 16 14:32:04 2022 +#> +#> Equations: +#> d_DMTA/dt = - k_DMTA * DMTA +#> d_M23/dt = + f_DMTA_to_M23 * k_DMTA * DMTA - k_M23 * M23 +#> d_M27/dt = + f_DMTA_to_M27 * k_DMTA * DMTA - k_M27 * M27 + k_M31 * M31 +#> d_M31/dt = + f_DMTA_to_M31 * k_DMTA * DMTA - k_M31 * M31 +#> +#> Data: +#> 563 observations of 4 variable(s) grouped in 6 datasets +#> +#> Model predictions using solution type deSolve +#> +#> Fitted in 927.963 s +#> Using 300, 100 iterations and 9 chains +#> +#> Variance model: Two-component variance function +#> +#> Mean of starting values for individual parameters: +#> DMTA_0 log_k_DMTA log_k_M23 log_k_M27 log_k_M31 f_DMTA_ilr_1 +#> 95.5662 -2.9048 -3.8130 -4.1600 -4.1486 0.1341 +#> f_DMTA_ilr_2 f_DMTA_ilr_3 +#> 0.1385 -1.6700 +#> +#> Fixed degradation parameter values: +#> None +#> +#> Results: +#> +#> Likelihood computed by importance sampling +#> AIC BIC logLik +#> 2276 2272 -1120 +#> +#> Optimised parameters: +#> est. lower upper +#> DMTA_0 88.5943 84.3961 92.7925 +#> log_k_DMTA -3.0466 -3.5609 -2.5322 +#> log_k_M23 -4.0684 -4.9340 -3.2029 +#> log_k_M27 -3.8628 -4.2627 -3.4628 +#> log_k_M31 -3.9803 -4.4804 -3.4801 +#> f_DMTA_ilr_1 0.1304 -0.2186 0.4795 +#> f_DMTA_ilr_2 0.1490 -0.2559 0.5540 +#> f_DMTA_ilr_3 -1.3970 -1.6976 -1.0964 +#> +#> Correlation: +#> DMTA_0 l__DMTA lg__M23 lg__M27 lg__M31 f_DMTA__1 f_DMTA__2 +#> log_k_DMTA 0.0309 +#> log_k_M23 -0.0231 -0.0031 +#> log_k_M27 -0.0381 -0.0048 0.0039 +#> log_k_M31 -0.0251 -0.0031 0.0021 0.0830 +#> f_DMTA_ilr_1 -0.0046 -0.0006 0.0417 -0.0437 0.0328 +#> f_DMTA_ilr_2 -0.0008 -0.0002 0.0214 -0.0270 -0.0909 -0.0361 +#> f_DMTA_ilr_3 -0.1832 -0.0135 0.0434 0.0804 0.0395 -0.0070 0.0059 +#> +#> Random effects: +#> est. lower upper +#> SD.DMTA_0 3.3651 -0.9655 7.6956 +#> SD.log_k_DMTA 0.6415 0.2774 1.0055 +#> SD.log_k_M23 1.0176 0.3809 1.6543 +#> SD.log_k_M27 0.4538 0.1522 0.7554 +#> SD.log_k_M31 0.5684 0.1905 0.9464 +#> SD.f_DMTA_ilr_1 0.4111 0.1524 0.6699 +#> SD.f_DMTA_ilr_2 0.4788 0.1808 0.7768 +#> SD.f_DMTA_ilr_3 0.3501 0.1316 0.5685 +#> +#> Variance model: +#> est. lower upper +#> a.1 0.9349 0.8395 1.0302 +#> b.1 0.1344 0.1176 0.1512 +#> +#> Backtransformed parameters: +#> est. lower upper +#> DMTA_0 88.59431 84.396147 92.79246 +#> k_DMTA 0.04752 0.028413 0.07948 +#> k_M23 0.01710 0.007198 0.04064 +#> k_M27 0.02101 0.014084 0.03134 +#> k_M31 0.01868 0.011329 0.03080 +#> f_DMTA_to_M23 0.14498 NA NA +#> f_DMTA_to_M27 0.12056 NA NA +#> f_DMTA_to_M31 0.11015 NA NA +#> +#> Resulting formation fractions: +#> ff +#> DMTA_M23 0.1450 +#> DMTA_M27 0.1206 +#> DMTA_M31 0.1101 +#> DMTA_sink 0.6243 +#> +#> Estimated disappearance times: +#> DT50 DT90 +#> DMTA 14.59 48.45 +#> M23 40.52 134.62 +#> M27 32.99 109.60 +#> M31 37.11 123.26 +# As the confidence interval for the random effects of DMTA_0 +# includes zero, we could try an alternative model without +# such random effects +# f_dmta_saem_tc_2 <- saem(dmta_sfo_sfo3p_tc, +# covariance.model = diag(c(0, rep(1, 7)))) +# saemix::plot(f_dmta_saem_tc_2$so, plot.type = "convergence") +# This does not perform better judged by AIC and BIC +saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so) +#> Error in saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so): object 'f_dmta_saem_tc_2' not found # }
    diff --git a/man/dimethenamid_2018.Rd b/man/dimethenamid_2018.Rd index 0d1265be..6c28ab7b 100644 --- a/man/dimethenamid_2018.Rd +++ b/man/dimethenamid_2018.Rd @@ -42,42 +42,43 @@ dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) dmta_ds[["Elliot 1"]] <- NULL dmta_ds[["Elliot 2"]] <- NULL \dontrun{ -dfop_sfo3_plus <- mkinmod( - DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), - M23 = mkinsub("SFO"), - M27 = mkinsub("SFO"), - M31 = mkinsub("SFO", "M27", sink = FALSE), - quiet = TRUE +# We don't use DFOP for the parent compound, as this gives numerical +# instabilities in the fits +sfo_sfo3p <- mkinmod( + DMTA = mkinsub("SFO", c("M23", "M27", "M31")), + M23 = mkinsub("SFO"), + M27 = mkinsub("SFO"), + M31 = mkinsub("SFO", "M27", sink = FALSE), + quiet = TRUE ) -f_dmta_mkin_tc <- mmkin( - list("DFOP-SFO3+" = dfop_sfo3_plus), - dmta_ds, quiet = TRUE, error_model = "tc") -nlmixr_model(f_dmta_mkin_tc) -# The focei fit takes about four minutes on my system -system.time( - f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", - control = nlmixr::foceiControl(print = 500)) -) -summary(f_dmta_nlmixr_focei) -plot(f_dmta_nlmixr_focei) -# Using saemix takes about 18 minutes -system.time( - f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) -) - -# nlmixr with est = "saem" is pretty fast with default iteration numbers, most -# of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end -# The likelihood calculated for the nlmixr fit is much lower than that found by saemix -# Also, the trace plot and the plot of the individual predictions is not -# convincing for the parent. It seems we are fitting an overparameterised -# model, so the result we get strongly depends on starting parameters and control settings. -system.time( - f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", - control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)) -) -traceplot(f_dmta_nlmixr_saem$nm) -summary(f_dmta_nlmixr_saem) -plot(f_dmta_nlmixr_saem) +dmta_sfo_sfo3p_tc <- mmkin(list("SFO-SFO3+" = sfo_sfo3p), + dmta_ds, error_model = "tc", quiet = TRUE) +print(dmta_sfo_sfo3p_tc) +# The default (test_log_parms = FALSE) gives an undue +# influence of ill-defined rate constants that have +# extremely small values: +plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = FALSE) +# If we disregards ill-defined rate constants, the results +# look more plausible, but the truth is likely to be in +# between these variants +plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = TRUE) +# Therefore we use nonlinear mixed-effects models +# f_dmta_nlme_tc <- nlme(dmta_sfo_sfo3p_tc) +# nlme reaches maxIter = 50 without convergence +f_dmta_saem_tc <- saem(dmta_sfo_sfo3p_tc) +# I am commenting out the convergence plot as rendering them +# with pkgdown fails (at least without further tweaks to the +# graphics device used) +#saemix::plot(f_dmta_saem_tc$so, plot.type = "convergence") +summary(f_dmta_saem_tc) +# As the confidence interval for the random effects of DMTA_0 +# includes zero, we could try an alternative model without +# such random effects +# f_dmta_saem_tc_2 <- saem(dmta_sfo_sfo3p_tc, +# covariance.model = diag(c(0, rep(1, 7)))) +# saemix::plot(f_dmta_saem_tc_2$so, plot.type = "convergence") +# This does not perform better judged by AIC and BIC +saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so) } } \keyword{datasets} diff --git a/vignettes/web_only/dimethenamid_2018.rmd b/vignettes/web_only/dimethenamid_2018.rmd index f86d776e..f34ca3ae 100644 --- a/vignettes/web_only/dimethenamid_2018.rmd +++ b/vignettes/web_only/dimethenamid_2018.rmd @@ -385,7 +385,7 @@ 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 in some for some of the soils. +cannot robustly be quantified for some of the soils. Regarding the regulatory use of these data, it is claimed that an improved characterisation of the mean parameter values across the population is -- cgit v1.2.1