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 --- docs/reference/dimethenamid_2018.html | 203 ++++++++++++++++++++++++---------- 1 file changed, 145 insertions(+), 58 deletions(-) (limited to 'docs/reference/dimethenamid_2018.html') 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 # } -- cgit v1.2.3