From b9be19af5e3085216d0cd5af439332f631fa8b92 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 15 Feb 2021 17:36:12 +0100 Subject: Fully rebuild docs, rerun tests and check --- docs/dev/reference/dimethenamid_2018.html | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'docs/dev/reference/dimethenamid_2018.html') diff --git a/docs/dev/reference/dimethenamid_2018.html b/docs/dev/reference/dimethenamid_2018.html index 21dea623..a06599df 100644 --- a/docs/dev/reference/dimethenamid_2018.html +++ b/docs/dev/reference/dimethenamid_2018.html @@ -77,7 +77,7 @@ constrained by data protection regulations." /> mkin - 0.9.50.4 + 1.0.3.9000 -- cgit v1.2.1 From 255430279d65bfe92093d48c9a586b062a38303d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 22 Jun 2021 15:15:02 +0200 Subject: Update development version of online docs --- docs/dev/reference/dimethenamid_2018.html | 288 +++++++++++++++++++++++++++++- 1 file changed, 286 insertions(+), 2 deletions(-) (limited to 'docs/dev/reference/dimethenamid_2018.html') diff --git a/docs/dev/reference/dimethenamid_2018.html b/docs/dev/reference/dimethenamid_2018.html index a06599df..e255765e 100644 --- a/docs/dev/reference/dimethenamid_2018.html +++ b/docs/dev/reference/dimethenamid_2018.html @@ -77,7 +77,7 @@ constrained by data protection regulations." /> mkin - 1.0.3.9000 + 1.0.5 @@ -203,7 +203,291 @@ specific pieces of information in the comments.

#> Elliot 2 0.75 33.37 23 #> Flaach 0.40 NA 20 #> BBA 2.2 0.40 NA 20 -#> BBA 2.3 0.40 NA 20 +#> BBA 2.3 0.40 NA 20
dmta_ds <- lapply(1:8, function(i) { + ds_i <- dimethenamid_2018$ds[[i]]$data + ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA" + ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i] + ds_i +}) +names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title) +dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]]) +dmta_ds[["Borstel 1"]] <- NULL +dmta_ds[["Borstel 2"]] <- NULL +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 +) +f_dmta_mkin_tc <- mmkin( + list("DFOP-SFO3+" = dfop_sfo3_plus), + dmta_ds, quiet = TRUE, error_model = "tc") +nlmixr_model(f_dmta_mkin_tc) +
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> function () +#> { +#> ini({ +#> DMTA_0 = 98.7697627680706 +#> eta.DMTA_0 ~ 2.35171765917765 +#> log_k_M23 = -3.92162409637283 +#> eta.log_k_M23 ~ 0.549278519419884 +#> log_k_M27 = -4.33774620773911 +#> eta.log_k_M27 ~ 0.864474956685295 +#> log_k_M31 = -4.24767627688461 +#> eta.log_k_M31 ~ 0.750297149164171 +#> log_k1 = -2.2341008812259 +#> eta.log_k1 ~ 0.902976221565793 +#> log_k2 = -3.7762779983269 +#> eta.log_k2 ~ 1.57684519529298 +#> g_qlogis = 0.450175725479389 +#> eta.g_qlogis ~ 3.0851335687675 +#> f_DMTA_tffm0_1_qlogis = -2.09240906629456 +#> eta.f_DMTA_tffm0_1_qlogis ~ 0.3 +#> f_DMTA_tffm0_2_qlogis = -2.18057573598794 +#> eta.f_DMTA_tffm0_2_qlogis ~ 0.3 +#> f_DMTA_tffm0_3_qlogis = -2.14267187609763 +#> eta.f_DMTA_tffm0_3_qlogis ~ 0.3 +#> sigma_low_DMTA = 0.697933852349996 +#> rsd_high_DMTA = 0.0257724286053519 +#> sigma_low_M23 = 0.697933852349996 +#> rsd_high_M23 = 0.0257724286053519 +#> sigma_low_M27 = 0.697933852349996 +#> rsd_high_M27 = 0.0257724286053519 +#> sigma_low_M31 = 0.697933852349996 +#> rsd_high_M31 = 0.0257724286053519 +#> }) +#> model({ +#> DMTA_0_model = DMTA_0 + eta.DMTA_0 +#> DMTA(0) = DMTA_0_model +#> k_M23 = exp(log_k_M23 + eta.log_k_M23) +#> k_M27 = exp(log_k_M27 + eta.log_k_M27) +#> k_M31 = exp(log_k_M31 + eta.log_k_M31) +#> k1 = exp(log_k1 + eta.log_k1) +#> k2 = exp(log_k2 + eta.log_k2) +#> g = expit(g_qlogis + eta.g_qlogis) +#> f_DMTA_tffm0_1 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis) +#> f_DMTA_tffm0_2 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis) +#> f_DMTA_tffm0_3 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis) +#> f_DMTA_to_M23 = f_DMTA_tffm0_1 +#> f_DMTA_to_M27 = f_DMTA_tffm0_2 * (1 - f_DMTA_tffm0_1) +#> f_DMTA_to_M31 = f_DMTA_tffm0_3 * (1 - f_DMTA_tffm0_2) * +#> (1 - f_DMTA_tffm0_1) +#> d/dt(DMTA) = -((k1 * g * exp(-k1 * time) + k2 * (1 - +#> g) * exp(-k2 * time))/(g * exp(-k1 * time) + (1 - +#> g) * exp(-k2 * time))) * DMTA +#> d/dt(M23) = +f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) + +#> k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) + +#> (1 - g) * exp(-k2 * time))) * DMTA - k_M23 * M23 +#> d/dt(M27) = +f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) + +#> k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) + +#> (1 - g) * exp(-k2 * time))) * DMTA - k_M27 * M27 + +#> k_M31 * M31 +#> d/dt(M31) = +f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) + +#> k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) + +#> (1 - g) * exp(-k2 * time))) * DMTA - k_M31 * M31 +#> DMTA ~ add(sigma_low_DMTA) + prop(rsd_high_DMTA) +#> M23 ~ add(sigma_low_M23) + prop(rsd_high_M23) +#> M27 ~ add(sigma_low_M27) + prop(rsd_high_M27) +#> M31 ~ add(sigma_low_M31) + prop(rsd_high_M31) +#> }) +#> } +#> <environment: 0x555559c2bd78>
f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", + control = nlmixr::foceiControl(print = 500)) +
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 +#>
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04 +#>
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01 +#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:08 +#>
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 +#>
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 +#>
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 +#>
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 +#>
#> → compiling inner model...
#>
#> ✔ done
#> → finding duplicate expressions in FD model...
#>
#> → optimizing duplicate expressions in FD model...
#>
#> → compiling EBE model...
#>
#> ✔ done
#> → compiling events FD model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> RxODE 1.1.0 using 8 threads (see ?getRxThreads) +#> no cache: create with `rxCreateCache()`
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation +#> F: Forward difference gradient approximation +#> C: Central difference gradient approximation +#> M: Mixed forward and central difference gradient approximation +#> Unscaled parameters for Omegas=chol(solve(omega)); +#> Diagonals are transformed, as specified by foceiControl(diagXform=) +#> |-----+---------------+-----------+-----------+-----------+-----------| +#> | #| Objective Fun | DMTA_0 | log_k_M23 | log_k_M27 | log_k_M31 | +#> |.....................| log_k1 | log_k2 | g_qlogis |f_DMTA_tffm0_1_qlogis | +#> |.....................|f_DMTA_tffm0_2_qlogis |f_DMTA_tffm0_3_qlogis | sigma_low | rsd_high | +#> |.....................| o1 | o2 | o3 | o4 | +#> |.....................| o5 | o6 | o7 | o8 | +#> |.....................| o9 | o10 |...........|...........| +#> calculating covariance matrix +#> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
summary(f_dmta_nlmixr_focei) +
#> nlmixr version used for fitting: 2.0.4 +#> mkin version used for pre-fitting: 1.0.5 +#> R version used for fitting: 4.1.0 +#> Date of fit: Thu Jun 17 14:04:58 2021 +#> Date of summary: Thu Jun 17 14:04:58 2021 +#> +#> Equations: +#> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * +#> time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) +#> * DMTA +#> d_M23/dt = + f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M23 * M23 +#> d_M27/dt = + f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M27 * M27 + k_M31 * M31 +#> d_M31/dt = + f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M31 * M31 +#> +#> Data: +#> 568 observations of 4 variable(s) grouped in 6 datasets +#> +#> Degradation model predictions using RxODE +#> +#> Fitted in 242.937 s +#> +#> Variance model: Two-component variance function +#> +#> Mean of starting values for individual parameters: +#> DMTA_0 log_k_M23 log_k_M27 log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 +#> 98.7698 -3.9216 -4.3377 -4.2477 0.1380 0.1393 +#> f_DMTA_ilr_3 log_k1 log_k2 g_qlogis +#> -1.7571 -2.2341 -3.7763 0.4502 +#> +#> Mean of starting values for error model parameters: +#> sigma_low rsd_high +#> 0.69793 0.02577 +#> +#> Fixed degradation parameter values: +#> None +#> +#> Results: +#> +#> Likelihood calculated by focei +#> AIC BIC logLik +#> 1936 2031 -945.9 +#> +#> Optimised parameters: +#> est. lower upper +#> DMTA_0 98.7698 98.7356 98.8039 +#> log_k_M23 -3.9216 -3.9235 -3.9197 +#> log_k_M27 -4.3377 -4.3398 -4.3357 +#> log_k_M31 -4.2477 -4.2497 -4.2457 +#> log_k1 -2.2341 -2.2353 -2.2329 +#> log_k2 -3.7763 -3.7781 -3.7744 +#> g_qlogis 0.4502 0.4496 0.4507 +#> f_DMTA_tffm0_1_qlogis -2.0924 -2.0936 -2.0912 +#> f_DMTA_tffm0_2_qlogis -2.1806 -2.1818 -2.1794 +#> f_DMTA_tffm0_3_qlogis -2.1427 -2.1439 -2.1415 +#> +#> Correlation: +#> DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs +#> log_k_M23 0 +#> log_k_M27 0 0 +#> log_k_M31 0 0 0 +#> log_k1 0 0 0 0 +#> log_k2 0 0 0 0 0 +#> g_qlogis 0 0 0 0 0 0 +#> f_DMTA_tffm0_1_qlogis 0 0 0 0 0 0 0 +#> f_DMTA_tffm0_2_qlogis 0 0 0 0 0 0 0 +#> f_DMTA_tffm0_3_qlogis 0 0 0 0 0 0 0 +#> f_DMTA_0_1 f_DMTA_0_2 +#> log_k_M23 +#> log_k_M27 +#> log_k_M31 +#> log_k1 +#> log_k2 +#> g_qlogis +#> f_DMTA_tffm0_1_qlogis +#> f_DMTA_tffm0_2_qlogis 0 +#> f_DMTA_tffm0_3_qlogis 0 0 +#> +#> Random effects (omega): +#> eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31 +#> eta.DMTA_0 2.352 0.0000 0.0000 0.0000 +#> eta.log_k_M23 0.000 0.5493 0.0000 0.0000 +#> eta.log_k_M27 0.000 0.0000 0.8645 0.0000 +#> eta.log_k_M31 0.000 0.0000 0.0000 0.7503 +#> eta.log_k1 0.000 0.0000 0.0000 0.0000 +#> eta.log_k2 0.000 0.0000 0.0000 0.0000 +#> eta.g_qlogis 0.000 0.0000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.0000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.0000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.0000 0.0000 0.0000 +#> eta.log_k1 eta.log_k2 eta.g_qlogis +#> eta.DMTA_0 0.000 0.000 0.000 +#> eta.log_k_M23 0.000 0.000 0.000 +#> eta.log_k_M27 0.000 0.000 0.000 +#> eta.log_k_M31 0.000 0.000 0.000 +#> eta.log_k1 0.903 0.000 0.000 +#> eta.log_k2 0.000 1.577 0.000 +#> eta.g_qlogis 0.000 0.000 3.085 +#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.000 0.000 +#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.000 0.000 +#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.000 0.000 +#> eta.f_DMTA_tffm0_1_qlogis eta.f_DMTA_tffm0_2_qlogis +#> eta.DMTA_0 0.0 0.0 +#> eta.log_k_M23 0.0 0.0 +#> eta.log_k_M27 0.0 0.0 +#> eta.log_k_M31 0.0 0.0 +#> eta.log_k1 0.0 0.0 +#> eta.log_k2 0.0 0.0 +#> eta.g_qlogis 0.0 0.0 +#> eta.f_DMTA_tffm0_1_qlogis 0.3 0.0 +#> eta.f_DMTA_tffm0_2_qlogis 0.0 0.3 +#> eta.f_DMTA_tffm0_3_qlogis 0.0 0.0 +#> eta.f_DMTA_tffm0_3_qlogis +#> eta.DMTA_0 0.0 +#> eta.log_k_M23 0.0 +#> eta.log_k_M27 0.0 +#> eta.log_k_M31 0.0 +#> eta.log_k1 0.0 +#> eta.log_k2 0.0 +#> eta.g_qlogis 0.0 +#> eta.f_DMTA_tffm0_1_qlogis 0.0 +#> eta.f_DMTA_tffm0_2_qlogis 0.0 +#> eta.f_DMTA_tffm0_3_qlogis 0.3 +#> +#> Variance model: +#> sigma_low rsd_high +#> 0.69793 0.02577 +#> +#> Backtransformed parameters: +#> est. lower upper +#> DMTA_0 98.76976 98.73563 98.80390 +#> k_M23 0.01981 0.01977 0.01985 +#> k_M27 0.01307 0.01304 0.01309 +#> k_M31 0.01430 0.01427 0.01433 +#> f_DMTA_to_M23 0.10984 NA NA +#> f_DMTA_to_M27 0.09036 NA NA +#> f_DMTA_to_M31 0.08399 NA NA +#> k1 0.10709 0.10696 0.10722 +#> k2 0.02291 0.02287 0.02295 +#> g 0.61068 0.61055 0.61081 +#> +#> Resulting formation fractions: +#> ff +#> DMTA_M23 0.10984 +#> DMTA_M27 0.09036 +#> DMTA_M31 0.08399 +#> DMTA_sink 0.71581 +#> +#> Estimated disappearance times: +#> DT50 DT90 DT50back DT50_k1 DT50_k2 +#> DMTA 10.66 59.78 18 6.473 30.26 +#> M23 34.99 116.24 NA NA NA +#> M27 53.05 176.23 NA NA NA +#> M31 48.48 161.05 NA NA NA
plot(f_dmta_nlmixr_focei) +
# saem has a problem with this model/data combination, maybe because of the +# overparameterised error model, to be investigated +#f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", +# control = saemControl(print = 500)) +#summary(f_dmta_nlmixr_saem) +#plot(f_dmta_nlmixr_saem) +# } +
@@ -168,7 +168,7 @@ constrained by data protection regulations.

Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018) Renewal Assessment Report Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour Rev. 2 - November 2017 -http://registerofquestions.efsa.europa.eu/roqFrontend/outputLoader?output=ON-5211

+https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716

Details

The R code used to create this data object is installed with this package @@ -295,8 +295,11 @@ specific pieces of information in the comments.

#> M31 ~ add(sigma_low_M31) + prop(rsd_high_M31) #> }) #> } -#> <environment: 0x555559c2bd78>
f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", - control = nlmixr::foceiControl(print = 500)) +#> <environment: 0x555559c00ce8>
# 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)) +)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 #>
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04 #>
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01 @@ -320,12 +323,13 @@ specific pieces of information in the comments.

#> |.....................| o5 | o6 | o7 | o8 | #> |.....................| o9 | o10 |...........|...........| #> calculating covariance matrix -#> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
summary(f_dmta_nlmixr_focei) +#> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
#> user system elapsed +#> 227.879 9.742 237.728
summary(f_dmta_nlmixr_focei)
#> nlmixr version used for fitting: 2.0.4 -#> mkin version used for pre-fitting: 1.0.5 +#> mkin version used for pre-fitting: 1.1.0 #> R version used for fitting: 4.1.0 -#> Date of fit: Thu Jun 17 14:04:58 2021 -#> Date of summary: Thu Jun 17 14:04:58 2021 +#> Date of fit: Tue Jul 27 16:02:33 2021 +#> Date of summary: Tue Jul 27 16:02:34 2021 #> #> Equations: #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -346,7 +350,7 @@ specific pieces of information in the comments.

#> #> Degradation model predictions using RxODE #> -#> Fitted in 242.937 s +#> Fitted in 237.547 s #> #> Variance model: Two-component variance function #> @@ -480,13 +484,194 @@ specific pieces of information in the comments.

#> M23 34.99 116.24 NA NA NA #> M27 53.05 176.23 NA NA NA #> M31 48.48 161.05 NA NA NA
plot(f_dmta_nlmixr_focei) -
# saem has a problem with this model/data combination, maybe because of the -# overparameterised error model, to be investigated -#f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", -# control = saemControl(print = 500)) -#summary(f_dmta_nlmixr_saem) -#plot(f_dmta_nlmixr_saem) -# } +
# Using saemix takes about 18 minutes +system.time( + f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) +) +
#> Running main SAEM algorithm +#> [1] "Tue Jul 27 16:02:34 2021" +#> .... +#> Minimisation finished +#> [1] "Tue Jul 27 16:21:39 2021"
#> user system elapsed +#> 1213.394 0.087 1213.578
+# 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)) +) +
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>
#> → generate SAEM model
#> ✔ done
#> 1: 98.3427 -3.5148 -3.3187 -3.7728 -2.1163 -2.8457 0.9482 -2.8064 -2.7412 -2.8745 2.7912 0.6805 0.8213 0.8055 0.8578 1.4980 2.9309 0.2850 0.2854 0.2850 4.0990 0.3821 3.5349 0.6537 5.4143 0.0002 4.5093 0.1905 +#> 500: 97.8277 -4.3506 -4.0318 -4.1520 -3.0553 -3.5843 1.1326 -2.0873 -2.0421 -2.0751 0.2960 1.2515 0.2531 0.3807 0.7928 0.8863 6.5211 0.1433 0.1082 0.3353 0.8960 0.0470 0.7501 0.0475 0.9527 0.0281 0.7321 0.0594
#> Calculating covariance matrix
#>
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#>
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#> user system elapsed +#> 818.782 3.808 154.926
traceplot(f_dmta_nlmixr_saem$nm) +
#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
summary(f_dmta_nlmixr_saem) +
#> nlmixr version used for fitting: 2.0.4 +#> mkin version used for pre-fitting: 1.1.0 +#> R version used for fitting: 4.1.0 +#> Date of fit: Tue Jul 27 16:25:23 2021 +#> Date of summary: Tue Jul 27 16:25:23 2021 +#> +#> Equations: +#> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * +#> time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) +#> * DMTA +#> d_M23/dt = + f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M23 * M23 +#> d_M27/dt = + f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M27 * M27 + k_M31 * M31 +#> d_M31/dt = + f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) +#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * +#> exp(-k2 * time))) * DMTA - k_M31 * M31 +#> +#> Data: +#> 568 observations of 4 variable(s) grouped in 6 datasets +#> +#> Degradation model predictions using RxODE +#> +#> Fitted in 154.632 s +#> +#> Variance model: Two-component variance function +#> +#> Mean of starting values for individual parameters: +#> DMTA_0 log_k_M23 log_k_M27 log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 +#> 98.7698 -3.9216 -4.3377 -4.2477 0.1380 0.1393 +#> f_DMTA_ilr_3 log_k1 log_k2 g_qlogis +#> -1.7571 -2.2341 -3.7763 0.4502 +#> +#> Mean of starting values for error model parameters: +#> sigma_low_DMTA rsd_high_DMTA sigma_low_M23 rsd_high_M23 sigma_low_M27 +#> 0.69793 0.02577 0.69793 0.02577 0.69793 +#> rsd_high_M27 sigma_low_M31 rsd_high_M31 +#> 0.02577 0.69793 0.02577 +#> +#> Fixed degradation parameter values: +#> None +#> +#> Results: +#> +#> Likelihood calculated by focei +#> AIC BIC logLik +#> 2036 2157 -989.8 +#> +#> Optimised parameters: +#> est. lower upper +#> DMTA_0 97.828 96.121 99.535 +#> log_k_M23 -4.351 -5.300 -3.401 +#> log_k_M27 -4.032 -4.470 -3.594 +#> log_k_M31 -4.152 -4.689 -3.615 +#> log_k1 -3.055 -3.785 -2.325 +#> log_k2 -3.584 -4.517 -2.651 +#> g_qlogis 1.133 -2.165 4.430 +#> f_DMTA_tffm0_1_qlogis -2.087 -2.407 -1.768 +#> f_DMTA_tffm0_2_qlogis -2.042 -2.336 -1.748 +#> f_DMTA_tffm0_3_qlogis -2.075 -2.557 -1.593 +#> +#> Correlation: +#> DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs +#> log_k_M23 -0.031 +#> log_k_M27 -0.050 0.004 +#> log_k_M31 -0.032 0.003 0.078 +#> log_k1 0.014 -0.002 -0.002 -0.001 +#> log_k2 0.059 0.006 -0.001 0.002 -0.037 +#> g_qlogis -0.077 0.005 0.009 0.004 0.035 -0.201 +#> f_DMTA_tffm0_1_qlogis -0.104 0.066 0.009 0.006 0.000 -0.011 0.014 +#> f_DMTA_tffm0_2_qlogis -0.120 0.013 0.081 -0.033 -0.002 -0.013 0.017 +#> f_DMTA_tffm0_3_qlogis -0.086 0.010 0.060 0.078 -0.002 -0.005 0.010 +#> f_DMTA_0_1 f_DMTA_0_2 +#> log_k_M23 +#> log_k_M27 +#> log_k_M31 +#> log_k1 +#> log_k2 +#> g_qlogis +#> f_DMTA_tffm0_1_qlogis +#> f_DMTA_tffm0_2_qlogis 0.026 +#> f_DMTA_tffm0_3_qlogis 0.019 0.002 +#> +#> Random effects (omega): +#> eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31 +#> eta.DMTA_0 0.296 0.000 0.0000 0.0000 +#> eta.log_k_M23 0.000 1.252 0.0000 0.0000 +#> eta.log_k_M27 0.000 0.000 0.2531 0.0000 +#> eta.log_k_M31 0.000 0.000 0.0000 0.3807 +#> eta.log_k1 0.000 0.000 0.0000 0.0000 +#> eta.log_k2 0.000 0.000 0.0000 0.0000 +#> eta.g_qlogis 0.000 0.000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.000 0.0000 0.0000 +#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.000 0.0000 0.0000 +#> eta.log_k1 eta.log_k2 eta.g_qlogis +#> eta.DMTA_0 0.0000 0.0000 0.000 +#> eta.log_k_M23 0.0000 0.0000 0.000 +#> eta.log_k_M27 0.0000 0.0000 0.000 +#> eta.log_k_M31 0.0000 0.0000 0.000 +#> eta.log_k1 0.7928 0.0000 0.000 +#> eta.log_k2 0.0000 0.8863 0.000 +#> eta.g_qlogis 0.0000 0.0000 6.521 +#> eta.f_DMTA_tffm0_1_qlogis 0.0000 0.0000 0.000 +#> eta.f_DMTA_tffm0_2_qlogis 0.0000 0.0000 0.000 +#> eta.f_DMTA_tffm0_3_qlogis 0.0000 0.0000 0.000 +#> eta.f_DMTA_tffm0_1_qlogis eta.f_DMTA_tffm0_2_qlogis +#> eta.DMTA_0 0.0000 0.0000 +#> eta.log_k_M23 0.0000 0.0000 +#> eta.log_k_M27 0.0000 0.0000 +#> eta.log_k_M31 0.0000 0.0000 +#> eta.log_k1 0.0000 0.0000 +#> eta.log_k2 0.0000 0.0000 +#> eta.g_qlogis 0.0000 0.0000 +#> eta.f_DMTA_tffm0_1_qlogis 0.1433 0.0000 +#> eta.f_DMTA_tffm0_2_qlogis 0.0000 0.1082 +#> eta.f_DMTA_tffm0_3_qlogis 0.0000 0.0000 +#> eta.f_DMTA_tffm0_3_qlogis +#> eta.DMTA_0 0.0000 +#> eta.log_k_M23 0.0000 +#> eta.log_k_M27 0.0000 +#> eta.log_k_M31 0.0000 +#> eta.log_k1 0.0000 +#> eta.log_k2 0.0000 +#> eta.g_qlogis 0.0000 +#> eta.f_DMTA_tffm0_1_qlogis 0.0000 +#> eta.f_DMTA_tffm0_2_qlogis 0.0000 +#> eta.f_DMTA_tffm0_3_qlogis 0.3353 +#> +#> Variance model: +#> sigma_low_DMTA rsd_high_DMTA sigma_low_M23 rsd_high_M23 sigma_low_M27 +#> 0.89603 0.04704 0.75015 0.04753 0.95265 +#> rsd_high_M27 sigma_low_M31 rsd_high_M31 +#> 0.02810 0.73212 0.05942 +#> +#> Backtransformed parameters: +#> est. lower upper +#> DMTA_0 97.82774 96.120503 99.53498 +#> k_M23 0.01290 0.004991 0.03334 +#> k_M27 0.01774 0.011451 0.02749 +#> k_M31 0.01573 0.009195 0.02692 +#> f_DMTA_to_M23 0.11033 NA NA +#> f_DMTA_to_M27 0.10218 NA NA +#> f_DMTA_to_M31 0.08784 NA NA +#> k1 0.04711 0.022707 0.09773 +#> k2 0.02775 0.010918 0.07056 +#> g 0.75632 0.102960 0.98823 +#> +#> Resulting formation fractions: +#> ff +#> DMTA_M23 0.11033 +#> DMTA_M27 0.10218 +#> DMTA_M31 0.08784 +#> DMTA_sink 0.69965 +#> +#> Estimated disappearance times: +#> DT50 DT90 DT50back DT50_k1 DT50_k2 +#> DMTA 16.59 57.44 17.29 14.71 24.97 +#> M23 53.74 178.51 NA NA NA +#> M27 39.07 129.78 NA NA NA +#> M31 44.06 146.36 NA NA NA
plot(f_dmta_nlmixr_saem) +
# }
@@ -295,7 +295,7 @@ specific pieces of information in the comments.

#> M31 ~ add(sigma_low_M31) + prop(rsd_high_M31) #> }) #> } -#> <environment: 0x555559c00ce8>
# The focei fit takes about four minutes on my system +#> <environment: 0x555559bfc940>
# 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)) @@ -303,7 +303,7 @@ specific pieces of information in the comments.

#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 #>
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04 #>
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01 -#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:08 +#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:09 #>
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 #>
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 #>
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 @@ -324,12 +324,12 @@ specific pieces of information in the comments.

#> |.....................| o9 | o10 |...........|...........| #> calculating covariance matrix #> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
#> user system elapsed -#> 227.879 9.742 237.728
summary(f_dmta_nlmixr_focei) +#> 227.223 8.444 235.624
summary(f_dmta_nlmixr_focei)
#> nlmixr version used for fitting: 2.0.4 -#> mkin version used for pre-fitting: 1.1.0 +#> mkin version used for pre-fitting: 1.0.5 #> R version used for fitting: 4.1.0 -#> Date of fit: Tue Jul 27 16:02:33 2021 -#> Date of summary: Tue Jul 27 16:02:34 2021 +#> Date of fit: Thu Jul 29 11:45:46 2021 +#> Date of summary: Thu Jul 29 11:45:46 2021 #> #> Equations: #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -350,7 +350,7 @@ specific pieces of information in the comments.

#> #> Degradation model predictions using RxODE #> -#> Fitted in 237.547 s +#> Fitted in 235.457 s #> #> Variance model: Two-component variance function #> @@ -489,11 +489,11 @@ specific pieces of information in the comments.

f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) )
#> Running main SAEM algorithm -#> [1] "Tue Jul 27 16:02:34 2021" +#> [1] "Thu Jul 29 11:45:47 2021" #> .... #> Minimisation finished -#> [1] "Tue Jul 27 16:21:39 2021"
#> user system elapsed -#> 1213.394 0.087 1213.578
+#> [1] "Thu Jul 29 12:04:25 2021"
#> user system elapsed +#> 1185.594 0.028 1185.687
# 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 @@ -506,13 +506,13 @@ specific pieces of information in the comments.

)
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>
#> → generate SAEM model
#> ✔ done
#> 1: 98.3427 -3.5148 -3.3187 -3.7728 -2.1163 -2.8457 0.9482 -2.8064 -2.7412 -2.8745 2.7912 0.6805 0.8213 0.8055 0.8578 1.4980 2.9309 0.2850 0.2854 0.2850 4.0990 0.3821 3.5349 0.6537 5.4143 0.0002 4.5093 0.1905 #> 500: 97.8277 -4.3506 -4.0318 -4.1520 -3.0553 -3.5843 1.1326 -2.0873 -2.0421 -2.0751 0.2960 1.2515 0.2531 0.3807 0.7928 0.8863 6.5211 0.1433 0.1082 0.3353 0.8960 0.0470 0.7501 0.0475 0.9527 0.0281 0.7321 0.0594
#> Calculating covariance matrix
#>
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#>
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#> user system elapsed -#> 818.782 3.808 154.926
traceplot(f_dmta_nlmixr_saem$nm) +#> 809.956 4.286 156.438
traceplot(f_dmta_nlmixr_saem$nm)
#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
summary(f_dmta_nlmixr_saem)
#> nlmixr version used for fitting: 2.0.4 -#> mkin version used for pre-fitting: 1.1.0 +#> mkin version used for pre-fitting: 1.0.5 #> R version used for fitting: 4.1.0 -#> Date of fit: Tue Jul 27 16:25:23 2021 -#> Date of summary: Tue Jul 27 16:25:23 2021 +#> Date of fit: Thu Jul 29 12:08:09 2021 +#> Date of summary: Thu Jul 29 12:08:09 2021 #> #> Equations: #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -533,7 +533,7 @@ specific pieces of information in the comments.

#> #> Degradation model predictions using RxODE #> -#> Fitted in 154.632 s +#> Fitted in 156.17 s #> #> Variance model: Two-component variance function #> -- cgit v1.2.1 From 51fab94230e926cec690dc455964bd797a97b7c7 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 4 Aug 2021 16:37:52 +0200 Subject: Improve AIC table in vignette --- docs/dev/reference/dimethenamid_2018.html | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) (limited to 'docs/dev/reference/dimethenamid_2018.html') diff --git a/docs/dev/reference/dimethenamid_2018.html b/docs/dev/reference/dimethenamid_2018.html index c893da63..a77cf0f4 100644 --- a/docs/dev/reference/dimethenamid_2018.html +++ b/docs/dev/reference/dimethenamid_2018.html @@ -295,7 +295,7 @@ specific pieces of information in the comments.

#> M31 ~ add(sigma_low_M31) + prop(rsd_high_M31) #> }) #> } -#> <environment: 0x555559bfc940>
# The focei fit takes about four minutes on my system +#> <environment: 0x555559ac3820>
# 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)) @@ -303,7 +303,7 @@ specific pieces of information in the comments.

#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 #>
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04 #>
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01 -#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:09 +#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:08 #>
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 #>
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 #>
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 @@ -324,12 +324,12 @@ specific pieces of information in the comments.

#> |.....................| o9 | o10 |...........|...........| #> calculating covariance matrix #> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
#> user system elapsed -#> 227.223 8.444 235.624
summary(f_dmta_nlmixr_focei) +#> 232.621 14.126 246.850
summary(f_dmta_nlmixr_focei)
#> nlmixr version used for fitting: 2.0.4 #> mkin version used for pre-fitting: 1.0.5 #> R version used for fitting: 4.1.0 -#> Date of fit: Thu Jul 29 11:45:46 2021 -#> Date of summary: Thu Jul 29 11:45:46 2021 +#> Date of fit: Wed Aug 4 15:53:54 2021 +#> Date of summary: Wed Aug 4 15:53:54 2021 #> #> Equations: #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -350,7 +350,7 @@ specific pieces of information in the comments.

#> #> Degradation model predictions using RxODE #> -#> Fitted in 235.457 s +#> Fitted in 246.669 s #> #> Variance model: Two-component variance function #> @@ -489,11 +489,11 @@ specific pieces of information in the comments.

f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) )
#> Running main SAEM algorithm -#> [1] "Thu Jul 29 11:45:47 2021" +#> [1] "Wed Aug 4 15:53:55 2021" #> .... #> Minimisation finished -#> [1] "Thu Jul 29 12:04:25 2021"
#> user system elapsed -#> 1185.594 0.028 1185.687
+#> [1] "Wed Aug 4 16:12:40 2021"
#> user system elapsed +#> 1192.021 0.064 1192.182
# 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 @@ -506,13 +506,13 @@ specific pieces of information in the comments.

)
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>
#> → generate SAEM model
#> ✔ done
#> 1: 98.3427 -3.5148 -3.3187 -3.7728 -2.1163 -2.8457 0.9482 -2.8064 -2.7412 -2.8745 2.7912 0.6805 0.8213 0.8055 0.8578 1.4980 2.9309 0.2850 0.2854 0.2850 4.0990 0.3821 3.5349 0.6537 5.4143 0.0002 4.5093 0.1905 #> 500: 97.8277 -4.3506 -4.0318 -4.1520 -3.0553 -3.5843 1.1326 -2.0873 -2.0421 -2.0751 0.2960 1.2515 0.2531 0.3807 0.7928 0.8863 6.5211 0.1433 0.1082 0.3353 0.8960 0.0470 0.7501 0.0475 0.9527 0.0281 0.7321 0.0594
#> Calculating covariance matrix
#>
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#>
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#> user system elapsed -#> 809.956 4.286 156.438
traceplot(f_dmta_nlmixr_saem$nm) +#> 813.299 3.736 151.935
traceplot(f_dmta_nlmixr_saem$nm)
#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
summary(f_dmta_nlmixr_saem)
#> nlmixr version used for fitting: 2.0.4 #> mkin version used for pre-fitting: 1.0.5 #> R version used for fitting: 4.1.0 -#> Date of fit: Thu Jul 29 12:08:09 2021 -#> Date of summary: Thu Jul 29 12:08:09 2021 +#> Date of fit: Wed Aug 4 16:16:18 2021 +#> Date of summary: Wed Aug 4 16:16:18 2021 #> #> Equations: #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -533,7 +533,7 @@ specific pieces of information in the comments.

#> #> Degradation model predictions using RxODE #> -#> Fitted in 156.17 s +#> Fitted in 151.67 s #> #> Variance model: Two-component variance function #> -- cgit v1.2.1 From c41381a961263c28d60976e68923157916c78b15 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 16 Sep 2021 15:31:13 +0200 Subject: Adapt and improve the dimethenamid vignette Adapt to the corrected data and unify control parameters for saemix and nlmixr with saem. Update docs --- docs/dev/reference/dimethenamid_2018.html | 320 +++++++++++++++--------------- 1 file changed, 157 insertions(+), 163 deletions(-) (limited to 'docs/dev/reference/dimethenamid_2018.html') diff --git a/docs/dev/reference/dimethenamid_2018.html b/docs/dev/reference/dimethenamid_2018.html index a77cf0f4..919e9363 100644 --- a/docs/dev/reference/dimethenamid_2018.html +++ b/docs/dev/reference/dimethenamid_2018.html @@ -77,7 +77,7 @@ constrained by data protection regulations." /> mkin - 1.0.5 + 1.1.0
@@ -162,7 +162,7 @@ constrained by data protection regulations.

Format

-

An mkindsg object grouping eight datasets with some meta information

+

An mkindsg object grouping seven datasets with some meta information

Source

Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018) @@ -177,42 +177,36 @@ specific pieces of information in the comments.

Examples

print(dimethenamid_2018)
-
#> <mkindsg> holding 8 mkinds objects
+
#> <mkindsg> holding 7 mkinds objects
 #> Title $title:  Aerobic soil degradation data on dimethenamid-P from the EU assessment in 2018 
 #> Occurrence of observed compounds $observed_n:
 #> DMTAP   M23   M27   M31  DMTA 
-#>     4     7     7     7     4 
+#>     3     7     7     7     4 
 #> Time normalisation factors $f_time_norm:
-#> [1] 1.0000000 0.9706477 0.9706477 1.2284784 1.2284784 0.6233856 0.7678922
-#> [8] 0.6733938
+#> [1] 1.0000000 0.9706477 1.2284784 1.2284784 0.6233856 0.7678922 0.6733938
 #> Meta information $meta:
-#>                       study  usda_soil_type study_moisture_ref_type
-#> Calke         Unsworth 2014      Sandy loam                     pF2
-#> Borstel 1 Staudenmaier 2013            Sand                     pF1
-#> Borstel 2 Staudenmaier 2009            Sand                     pF1
-#> Elliot 1         Wendt 1997       Clay loam                   pF2.5
-#> Elliot 2         Wendt 1997       Clay loam                   pF2.5
-#> Flaach           König 1996 Sandy clay loam                     pF1
-#> BBA 2.2          König 1995      Loamy sand                     pF1
-#> BBA 2.3          König 1995      Sandy loam                     pF1
-#>           rel_moisture study_ref_moisture temperature
-#> Calke             1.00                 NA          20
-#> Borstel 1         0.50              23.00          20
-#> Borstel 2         0.50              23.00          20
-#> Elliot 1          0.75              33.37          23
-#> Elliot 2          0.75              33.37          23
-#> Flaach            0.40                 NA          20
-#> BBA 2.2           0.40                 NA          20
-#> BBA 2.3           0.40                 NA          20
dmta_ds <- lapply(1:8, function(i) {
+#>                      study  usda_soil_type study_moisture_ref_type rel_moisture
+#> Calke        Unsworth 2014      Sandy loam                     pF2         1.00
+#> Borstel  Staudenmaier 2009            Sand                     pF1         0.50
+#> Elliot 1        Wendt 1997       Clay loam                   pF2.5         0.75
+#> Elliot 2        Wendt 1997       Clay loam                   pF2.5         0.75
+#> Flaach          König 1996 Sandy clay loam                     pF1         0.40
+#> BBA 2.2         König 1995      Loamy sand                     pF1         0.40
+#> BBA 2.3         König 1995      Sandy loam                     pF1         0.40
+#>          study_ref_moisture temperature
+#> Calke                    NA          20
+#> Borstel               23.00          20
+#> Elliot 1              33.37          23
+#> Elliot 2              33.37          23
+#> Flaach                   NA          20
+#> BBA 2.2                  NA          20
+#> BBA 2.3                  NA          20
dmta_ds <- lapply(1:7, function(i) {
   ds_i <- dimethenamid_2018$ds[[i]]$data
   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
   ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
   ds_i
 })
 names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
-dmta_ds[["Borstel 1"]] <- NULL
-dmta_ds[["Borstel 2"]] <- NULL
 dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
 dmta_ds[["Elliot 1"]] <- NULL
 dmta_ds[["Elliot 2"]] <- NULL
@@ -231,33 +225,33 @@ specific pieces of information in the comments.

 
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> function () 
 #> {
 #>     ini({
-#>         DMTA_0 = 98.7697627680706
-#>         eta.DMTA_0 ~ 2.35171765917765
+#>         DMTA_0 = 98.7132391714013
+#>         eta.DMTA_0 ~ 2.32692496033921
 #>         log_k_M23 = -3.92162409637283
 #>         eta.log_k_M23 ~ 0.549278519419884
-#>         log_k_M27 = -4.33774620773911
-#>         eta.log_k_M27 ~ 0.864474956685295
-#>         log_k_M31 = -4.24767627688461
-#>         eta.log_k_M31 ~ 0.750297149164171
-#>         log_k1 = -2.2341008812259
-#>         eta.log_k1 ~ 0.902976221565793
-#>         log_k2 = -3.7762779983269
-#>         eta.log_k2 ~ 1.57684519529298
-#>         g_qlogis = 0.450175725479389
-#>         eta.g_qlogis ~ 3.0851335687675
-#>         f_DMTA_tffm0_1_qlogis = -2.09240906629456
+#>         log_k_M27 = -4.33057580082049
+#>         eta.log_k_M27 ~ 0.855184233768426
+#>         log_k_M31 = -4.24415516780733
+#>         eta.log_k_M31 ~ 0.745746058085877
+#>         log_k1 = -2.23515804885306
+#>         eta.log_k1 ~ 0.901033446532357
+#>         log_k2 = -3.77581484944379
+#>         eta.log_k2 ~ 1.57682329638124
+#>         g_qlogis = 0.436302910942805
+#>         eta.g_qlogis ~ 3.10190528862808
+#>         f_DMTA_tffm0_1_qlogis = -2.0914852208395
 #>         eta.f_DMTA_tffm0_1_qlogis ~ 0.3
-#>         f_DMTA_tffm0_2_qlogis = -2.18057573598794
+#>         f_DMTA_tffm0_2_qlogis = -2.17879574608926
 #>         eta.f_DMTA_tffm0_2_qlogis ~ 0.3
-#>         f_DMTA_tffm0_3_qlogis = -2.14267187609763
+#>         f_DMTA_tffm0_3_qlogis = -2.14036526460782
 #>         eta.f_DMTA_tffm0_3_qlogis ~ 0.3
-#>         sigma_low_DMTA = 0.697933852349996
+#>         sigma_low_DMTA = 0.700117227383809
 #>         rsd_high_DMTA = 0.0257724286053519
-#>         sigma_low_M23 = 0.697933852349996
+#>         sigma_low_M23 = 0.700117227383809
 #>         rsd_high_M23 = 0.0257724286053519
-#>         sigma_low_M27 = 0.697933852349996
+#>         sigma_low_M27 = 0.700117227383809
 #>         rsd_high_M27 = 0.0257724286053519
-#>         sigma_low_M31 = 0.697933852349996
+#>         sigma_low_M31 = 0.700117227383809
 #>         rsd_high_M31 = 0.0257724286053519
 #>     })
 #>     model({
@@ -295,7 +289,7 @@ specific pieces of information in the comments.

 #>         M31 ~ add(sigma_low_M31) + prop(rsd_high_M31)
 #>     })
 #> }
-#> <environment: 0x555559ac3820>
# The focei fit takes about four minutes on my system
+#> <environment: 0x555559e97ac0>
# 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))
@@ -308,7 +302,7 @@ specific pieces of information in the comments.

 #> 
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 
 #> 
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
 #> 
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
-#> 
#> → compiling inner model...
#>  
#> ✔ done
#> → finding duplicate expressions in FD model...
#> 
#> → optimizing duplicate expressions in FD model...
#> 
#> → compiling EBE model...
#>  
#> ✔ done
#> → compiling events FD model...
#>  
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> RxODE 1.1.0 using 8 threads (see ?getRxThreads)
+#> 
#> → compiling inner model...
#>  
#> ✔ done
#> → finding duplicate expressions in FD model...
#> 
#> → optimizing duplicate expressions in FD model...
#> 
#> → compiling EBE model...
#>  
#> ✔ done
#> → compiling events FD model...
#>  
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> RxODE 1.1.1 using 8 threads (see ?getRxThreads)
 #>   no cache: create with `rxCreateCache()`
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
 #> F: Forward difference gradient approximation
 #> C: Central difference gradient approximation
@@ -324,12 +318,12 @@ specific pieces of information in the comments.

 #> |.....................|        o9 |       o10 |...........|...........|
 #> calculating covariance matrix
 #> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
#>    user  system elapsed 
-#> 232.621  14.126 246.850 
summary(f_dmta_nlmixr_focei)
-
#> nlmixr version used for fitting:    2.0.4 
-#> mkin version used for pre-fitting:  1.0.5 
-#> R version used for fitting:         4.1.0 
-#> Date of fit:     Wed Aug  4 15:53:54 2021 
-#> Date of summary: Wed Aug  4 15:53:54 2021 
+#> 230.015   8.962 238.957 
summary(f_dmta_nlmixr_focei)
+
#> nlmixr version used for fitting:    2.0.5 
+#> mkin version used for pre-fitting:  1.1.0 
+#> R version used for fitting:         4.1.1 
+#> Date of fit:     Thu Sep 16 14:06:55 2021 
+#> Date of summary: Thu Sep 16 14:06:55 2021 
 #> 
 #> Equations:
 #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -346,23 +340,23 @@ specific pieces of information in the comments.

 #>            exp(-k2 * time))) * DMTA - k_M31 * M31
 #> 
 #> Data:
-#> 568 observations of 4 variable(s) grouped in 6 datasets
+#> 563 observations of 4 variable(s) grouped in 6 datasets
 #> 
 #> Degradation model predictions using RxODE
 #> 
-#> Fitted in 246.669 s
+#> Fitted in 238.792 s
 #> 
 #> Variance model: Two-component variance function 
 #> 
 #> Mean of starting values for individual parameters:
 #>       DMTA_0    log_k_M23    log_k_M27    log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 
-#>      98.7698      -3.9216      -4.3377      -4.2477       0.1380       0.1393 
+#>      98.7132      -3.9216      -4.3306      -4.2442       0.1376       0.1388 
 #> f_DMTA_ilr_3       log_k1       log_k2     g_qlogis 
-#>      -1.7571      -2.2341      -3.7763       0.4502 
+#>      -1.7554      -2.2352      -3.7758       0.4363 
 #> 
 #> Mean of starting values for error model parameters:
 #> sigma_low  rsd_high 
-#>   0.69793   0.02577 
+#>   0.70012   0.02577 
 #> 
 #> Fixed degradation parameter values:
 #> None
@@ -371,20 +365,20 @@ specific pieces of information in the comments.

 #> 
 #> Likelihood calculated by focei  
 #>    AIC  BIC logLik
-#>   1936 2031 -945.9
+#>   1918 2014 -937.2
 #> 
 #> Optimised parameters:
 #>                          est.   lower   upper
-#> DMTA_0                98.7698 98.7356 98.8039
-#> log_k_M23             -3.9216 -3.9235 -3.9197
-#> log_k_M27             -4.3377 -4.3398 -4.3357
-#> log_k_M31             -4.2477 -4.2497 -4.2457
-#> log_k1                -2.2341 -2.2353 -2.2329
-#> log_k2                -3.7763 -3.7781 -3.7744
-#> g_qlogis               0.4502  0.4496  0.4507
-#> f_DMTA_tffm0_1_qlogis -2.0924 -2.0936 -2.0912
-#> f_DMTA_tffm0_2_qlogis -2.1806 -2.1818 -2.1794
-#> f_DMTA_tffm0_3_qlogis -2.1427 -2.1439 -2.1415
+#> DMTA_0                98.7132 98.6801 98.7464
+#> log_k_M23             -3.9216 -3.9235 -3.9198
+#> log_k_M27             -4.3306 -4.3326 -4.3286
+#> log_k_M31             -4.2442 -4.2461 -4.2422
+#> log_k1                -2.2352 -2.2364 -2.2340
+#> log_k2                -3.7758 -3.7776 -3.7740
+#> g_qlogis               0.4363  0.4358  0.4368
+#> f_DMTA_tffm0_1_qlogis -2.0915 -2.0926 -2.0903
+#> f_DMTA_tffm0_2_qlogis -2.1788 -2.1800 -2.1776
+#> f_DMTA_tffm0_3_qlogis -2.1404 -2.1415 -2.1392
 #> 
 #> Correlation: 
 #>                       DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs
@@ -410,10 +404,10 @@ specific pieces of information in the comments.

 #> 
 #> Random effects (omega):
 #>                           eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31
-#> eta.DMTA_0                     2.352        0.0000        0.0000        0.0000
+#> eta.DMTA_0                     2.327        0.0000        0.0000        0.0000
 #> eta.log_k_M23                  0.000        0.5493        0.0000        0.0000
-#> eta.log_k_M27                  0.000        0.0000        0.8645        0.0000
-#> eta.log_k_M31                  0.000        0.0000        0.0000        0.7503
+#> eta.log_k_M27                  0.000        0.0000        0.8552        0.0000
+#> eta.log_k_M31                  0.000        0.0000        0.0000        0.7457
 #> eta.log_k1                     0.000        0.0000        0.0000        0.0000
 #> eta.log_k2                     0.000        0.0000        0.0000        0.0000
 #> eta.g_qlogis                   0.000        0.0000        0.0000        0.0000
@@ -425,9 +419,9 @@ specific pieces of information in the comments.

 #> eta.log_k_M23                  0.000      0.000        0.000
 #> eta.log_k_M27                  0.000      0.000        0.000
 #> eta.log_k_M31                  0.000      0.000        0.000
-#> eta.log_k1                     0.903      0.000        0.000
+#> eta.log_k1                     0.901      0.000        0.000
 #> eta.log_k2                     0.000      1.577        0.000
-#> eta.g_qlogis                   0.000      0.000        3.085
+#> eta.g_qlogis                   0.000      0.000        3.102
 #> eta.f_DMTA_tffm0_1_qlogis      0.000      0.000        0.000
 #> eta.f_DMTA_tffm0_2_qlogis      0.000      0.000        0.000
 #> eta.f_DMTA_tffm0_3_qlogis      0.000      0.000        0.000
@@ -456,44 +450,44 @@ specific pieces of information in the comments.

 #> 
 #> Variance model:
 #> sigma_low  rsd_high 
-#>   0.69793   0.02577 
+#>   0.70012   0.02577 
 #> 
 #> Backtransformed parameters:
 #>                   est.    lower    upper
-#> DMTA_0        98.76976 98.73563 98.80390
+#> DMTA_0        98.71324 98.68012 98.74636
 #> k_M23          0.01981  0.01977  0.01985
-#> k_M27          0.01307  0.01304  0.01309
-#> k_M31          0.01430  0.01427  0.01433
-#> f_DMTA_to_M23  0.10984       NA       NA
-#> f_DMTA_to_M27  0.09036       NA       NA
-#> f_DMTA_to_M31  0.08399       NA       NA
-#> k1             0.10709  0.10696  0.10722
-#> k2             0.02291  0.02287  0.02295
-#> g              0.61068  0.61055  0.61081
+#> k_M27          0.01316  0.01313  0.01319
+#> k_M31          0.01435  0.01432  0.01438
+#> f_DMTA_to_M23  0.10993       NA       NA
+#> f_DMTA_to_M27  0.09049       NA       NA
+#> f_DMTA_to_M31  0.08414       NA       NA
+#> k1             0.10698  0.10685  0.10710
+#> k2             0.02292  0.02288  0.02296
+#> g              0.60738  0.60725  0.60751
 #> 
 #> Resulting formation fractions:
 #>                ff
-#> DMTA_M23  0.10984
-#> DMTA_M27  0.09036
-#> DMTA_M31  0.08399
-#> DMTA_sink 0.71581
+#> DMTA_M23  0.10993
+#> DMTA_M27  0.09049
+#> DMTA_M31  0.08414
+#> DMTA_sink 0.71543
 #> 
 #> Estimated disappearance times:
-#>       DT50   DT90 DT50back DT50_k1 DT50_k2
-#> DMTA 10.66  59.78       18   6.473   30.26
-#> M23  34.99 116.24       NA      NA      NA
-#> M27  53.05 176.23       NA      NA      NA
-#> M31  48.48 161.05       NA      NA      NA
plot(f_dmta_nlmixr_focei)
+#>       DT50  DT90 DT50back DT50_k1 DT50_k2
+#> DMTA 10.72  60.1    18.09    6.48   30.24
+#> M23  34.99 116.2       NA      NA      NA
+#> M27  52.67 175.0       NA      NA      NA
+#> M31  48.31 160.5       NA      NA      NA
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)
 )
 
#> Running main SAEM algorithm
-#> [1] "Wed Aug  4 15:53:55 2021"
+#> [1] "Thu Sep 16 14:06:56 2021"
 #> ....
 #>     Minimisation finished
-#> [1] "Wed Aug  4 16:12:40 2021"
#>     user   system  elapsed 
-#> 1192.021    0.064 1192.182 

+#> [1] "Thu Sep 16 14:25:28 2021"
#>     user   system  elapsed 
+#> 1176.278    0.021 1176.388 

 # 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
@@ -504,15 +498,15 @@ specific pieces of information in the comments.

   f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
     control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9))
 )
-
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>  
#> → generate SAEM model
#> ✔ done
#> 1:    98.3427   -3.5148   -3.3187   -3.7728   -2.1163   -2.8457    0.9482   -2.8064   -2.7412   -2.8745    2.7912    0.6805    0.8213    0.8055    0.8578    1.4980    2.9309    0.2850    0.2854    0.2850    4.0990    0.3821    3.5349    0.6537    5.4143    0.0002    4.5093    0.1905
-#> 500:    97.8277   -4.3506   -4.0318   -4.1520   -3.0553   -3.5843    1.1326   -2.0873   -2.0421   -2.0751    0.2960    1.2515    0.2531    0.3807    0.7928    0.8863    6.5211    0.1433    0.1082    0.3353    0.8960    0.0470    0.7501    0.0475    0.9527    0.0281    0.7321    0.0594
#> Calculating covariance matrix
#> 
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#> 
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>  
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#>    user  system elapsed 
-#> 813.299   3.736 151.935 
traceplot(f_dmta_nlmixr_saem$nm)
+
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>  
#> → generate SAEM model
#> ✔ done
#> 1:    98.3400   -3.5096   -3.3392   -3.7596   -2.2055   -2.7755    1.0281   -2.7872   -2.7223   -2.8341    2.6422    0.7027    0.8124    0.7085    0.8560    1.4980    3.2777    0.3063    0.2850    0.2850    4.1120    0.3716    4.4582    0.3994    4.4820    0.4025    3.7803    0.5780
+#> 500:    97.8212   -4.4030   -4.0872   -4.1289   -2.8278   -4.3505    2.6614   -2.1252   -2.1308   -2.0749    2.9463    1.2933    0.2802    0.3467    0.4814    0.7877    3.0743    0.1508    0.1523    0.3155    0.9557    0.0333    0.4787    0.1073    0.6826    0.0707    0.7849    0.0356
#> Calculating covariance matrix
#> 
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#> 
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>  
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#>    user  system elapsed 
+#> 800.784   3.715 149.687 
traceplot(f_dmta_nlmixr_saem$nm)
 
#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
summary(f_dmta_nlmixr_saem)
-
#> nlmixr version used for fitting:    2.0.4 
-#> mkin version used for pre-fitting:  1.0.5 
-#> R version used for fitting:         4.1.0 
-#> Date of fit:     Wed Aug  4 16:16:18 2021 
-#> Date of summary: Wed Aug  4 16:16:18 2021 
+
#> nlmixr version used for fitting:    2.0.5 
+#> mkin version used for pre-fitting:  1.1.0 
+#> R version used for fitting:         4.1.1 
+#> Date of fit:     Thu Sep 16 14:29:02 2021 
+#> Date of summary: Thu Sep 16 14:29:02 2021 
 #> 
 #> Equations:
 #> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -529,25 +523,25 @@ specific pieces of information in the comments.

 #>            exp(-k2 * time))) * DMTA - k_M31 * M31
 #> 
 #> Data:
-#> 568 observations of 4 variable(s) grouped in 6 datasets
+#> 563 observations of 4 variable(s) grouped in 6 datasets
 #> 
 #> Degradation model predictions using RxODE
 #> 
-#> Fitted in 151.67 s
+#> Fitted in 149.421 s
 #> 
 #> Variance model: Two-component variance function 
 #> 
 #> Mean of starting values for individual parameters:
 #>       DMTA_0    log_k_M23    log_k_M27    log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 
-#>      98.7698      -3.9216      -4.3377      -4.2477       0.1380       0.1393 
+#>      98.7132      -3.9216      -4.3306      -4.2442       0.1376       0.1388 
 #> f_DMTA_ilr_3       log_k1       log_k2     g_qlogis 
-#>      -1.7571      -2.2341      -3.7763       0.4502 
+#>      -1.7554      -2.2352      -3.7758       0.4363 
 #> 
 #> Mean of starting values for error model parameters:
 #> sigma_low_DMTA  rsd_high_DMTA  sigma_low_M23   rsd_high_M23  sigma_low_M27 
-#>        0.69793        0.02577        0.69793        0.02577        0.69793 
+#>        0.70012        0.02577        0.70012        0.02577        0.70012 
 #>   rsd_high_M27  sigma_low_M31   rsd_high_M31 
-#>        0.02577        0.69793        0.02577 
+#>        0.02577        0.70012        0.02577 
 #> 
 #> Fixed degradation parameter values:
 #> None
@@ -556,32 +550,32 @@ specific pieces of information in the comments.

 #> 
 #> Likelihood calculated by focei  
 #>    AIC  BIC logLik
-#>   2036 2157 -989.8
+#>   1953 2074 -948.3
 #> 
 #> Optimised parameters:
 #>                         est.  lower  upper
-#> DMTA_0                97.828 96.121 99.535
-#> log_k_M23             -4.351 -5.300 -3.401
-#> log_k_M27             -4.032 -4.470 -3.594
-#> log_k_M31             -4.152 -4.689 -3.615
-#> log_k1                -3.055 -3.785 -2.325
-#> log_k2                -3.584 -4.517 -2.651
-#> g_qlogis               1.133 -2.165  4.430
-#> f_DMTA_tffm0_1_qlogis -2.087 -2.407 -1.768
-#> f_DMTA_tffm0_2_qlogis -2.042 -2.336 -1.748
-#> f_DMTA_tffm0_3_qlogis -2.075 -2.557 -1.593
+#> DMTA_0                97.821 95.862 99.780
+#> log_k_M23             -4.403 -5.376 -3.430
+#> log_k_M27             -4.087 -4.545 -3.629
+#> log_k_M31             -4.129 -4.639 -3.618
+#> log_k1                -2.828 -3.389 -2.266
+#> log_k2                -4.351 -5.472 -3.229
+#> g_qlogis               2.661  0.824  4.499
+#> f_DMTA_tffm0_1_qlogis -2.125 -2.449 -1.801
+#> f_DMTA_tffm0_2_qlogis -2.131 -2.468 -1.794
+#> f_DMTA_tffm0_3_qlogis -2.075 -2.540 -1.610
 #> 
 #> Correlation: 
 #>                       DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs
-#> log_k_M23             -0.031                                          
-#> log_k_M27             -0.050  0.004                                   
-#> log_k_M31             -0.032  0.003  0.078                            
-#> log_k1                 0.014 -0.002 -0.002 -0.001                     
-#> log_k2                 0.059  0.006 -0.001  0.002 -0.037              
-#> g_qlogis              -0.077  0.005  0.009  0.004  0.035 -0.201       
-#> f_DMTA_tffm0_1_qlogis -0.104  0.066  0.009  0.006  0.000 -0.011  0.014
-#> f_DMTA_tffm0_2_qlogis -0.120  0.013  0.081 -0.033 -0.002 -0.013  0.017
-#> f_DMTA_tffm0_3_qlogis -0.086  0.010  0.060  0.078 -0.002 -0.005  0.010
+#> log_k_M23             -0.019                                          
+#> log_k_M27             -0.028  0.004                                   
+#> log_k_M31             -0.019  0.003  0.075                            
+#> log_k1                 0.038 -0.004 -0.006 -0.003                     
+#> log_k2                 0.046  0.011  0.008  0.009  0.068              
+#> g_qlogis              -0.067  0.004  0.006  0.001 -0.076 -0.409       
+#> f_DMTA_tffm0_1_qlogis -0.062  0.055  0.006  0.004 -0.008 -0.004  0.012
+#> f_DMTA_tffm0_2_qlogis -0.062  0.010  0.058 -0.034 -0.008 -0.007  0.014
+#> f_DMTA_tffm0_3_qlogis -0.052  0.009  0.056  0.071 -0.006 -0.001  0.008
 #>                       f_DMTA_0_1 f_DMTA_0_2
 #> log_k_M23                                  
 #> log_k_M27                                  
@@ -590,15 +584,15 @@ specific pieces of information in the comments.

 #> log_k2                                     
 #> g_qlogis                                   
 #> f_DMTA_tffm0_1_qlogis                      
-#> f_DMTA_tffm0_2_qlogis  0.026               
-#> f_DMTA_tffm0_3_qlogis  0.019      0.002    
+#> f_DMTA_tffm0_2_qlogis  0.017               
+#> f_DMTA_tffm0_3_qlogis  0.014     -0.005    
 #> 
 #> Random effects (omega):
 #>                           eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31
-#> eta.DMTA_0                     0.296         0.000        0.0000        0.0000
-#> eta.log_k_M23                  0.000         1.252        0.0000        0.0000
-#> eta.log_k_M27                  0.000         0.000        0.2531        0.0000
-#> eta.log_k_M31                  0.000         0.000        0.0000        0.3807
+#> eta.DMTA_0                     2.946         0.000        0.0000        0.0000
+#> eta.log_k_M23                  0.000         1.293        0.0000        0.0000
+#> eta.log_k_M27                  0.000         0.000        0.2802        0.0000
+#> eta.log_k_M31                  0.000         0.000        0.0000        0.3467
 #> eta.log_k1                     0.000         0.000        0.0000        0.0000
 #> eta.log_k2                     0.000         0.000        0.0000        0.0000
 #> eta.g_qlogis                   0.000         0.000        0.0000        0.0000
@@ -610,9 +604,9 @@ specific pieces of information in the comments.

 #> eta.log_k_M23                 0.0000     0.0000        0.000
 #> eta.log_k_M27                 0.0000     0.0000        0.000
 #> eta.log_k_M31                 0.0000     0.0000        0.000
-#> eta.log_k1                    0.7928     0.0000        0.000
-#> eta.log_k2                    0.0000     0.8863        0.000
-#> eta.g_qlogis                  0.0000     0.0000        6.521
+#> eta.log_k1                    0.4814     0.0000        0.000
+#> eta.log_k2                    0.0000     0.7877        0.000
+#> eta.g_qlogis                  0.0000     0.0000        3.074
 #> eta.f_DMTA_tffm0_1_qlogis     0.0000     0.0000        0.000
 #> eta.f_DMTA_tffm0_2_qlogis     0.0000     0.0000        0.000
 #> eta.f_DMTA_tffm0_3_qlogis     0.0000     0.0000        0.000
@@ -624,8 +618,8 @@ specific pieces of information in the comments.

 #> eta.log_k1                                   0.0000                    0.0000
 #> eta.log_k2                                   0.0000                    0.0000
 #> eta.g_qlogis                                 0.0000                    0.0000
-#> eta.f_DMTA_tffm0_1_qlogis                    0.1433                    0.0000
-#> eta.f_DMTA_tffm0_2_qlogis                    0.0000                    0.1082
+#> eta.f_DMTA_tffm0_1_qlogis                    0.1508                    0.0000
+#> eta.f_DMTA_tffm0_2_qlogis                    0.0000                    0.1523
 #> eta.f_DMTA_tffm0_3_qlogis                    0.0000                    0.0000
 #>                           eta.f_DMTA_tffm0_3_qlogis
 #> eta.DMTA_0                                   0.0000
@@ -637,40 +631,40 @@ specific pieces of information in the comments.

 #> eta.g_qlogis                                 0.0000
 #> eta.f_DMTA_tffm0_1_qlogis                    0.0000
 #> eta.f_DMTA_tffm0_2_qlogis                    0.0000
-#> eta.f_DMTA_tffm0_3_qlogis                    0.3353
+#> eta.f_DMTA_tffm0_3_qlogis                    0.3155
 #> 
 #> Variance model:
 #> sigma_low_DMTA  rsd_high_DMTA  sigma_low_M23   rsd_high_M23  sigma_low_M27 
-#>        0.89603        0.04704        0.75015        0.04753        0.95265 
+#>        0.95572        0.03325        0.47871        0.10733        0.68264 
 #>   rsd_high_M27  sigma_low_M31   rsd_high_M31 
-#>        0.02810        0.73212        0.05942 
+#>        0.07072        0.78486        0.03557 
 #> 
 #> Backtransformed parameters:
 #>                   est.     lower    upper
-#> DMTA_0        97.82774 96.120503 99.53498
-#> k_M23          0.01290  0.004991  0.03334
-#> k_M27          0.01774  0.011451  0.02749
-#> k_M31          0.01573  0.009195  0.02692
-#> f_DMTA_to_M23  0.11033        NA       NA
-#> f_DMTA_to_M27  0.10218        NA       NA
-#> f_DMTA_to_M31  0.08784        NA       NA
-#> k1             0.04711  0.022707  0.09773
-#> k2             0.02775  0.010918  0.07056
-#> g              0.75632  0.102960  0.98823
+#> DMTA_0        97.82122 95.862233 99.78020
+#> k_M23          0.01224  0.004625  0.03239
+#> k_M27          0.01679  0.010615  0.02654
+#> k_M31          0.01610  0.009664  0.02683
+#> f_DMTA_to_M23  0.10668        NA       NA
+#> f_DMTA_to_M27  0.09481        NA       NA
+#> f_DMTA_to_M31  0.08908        NA       NA
+#> k1             0.05914  0.033731  0.10370
+#> k2             0.01290  0.004204  0.03958
+#> g              0.93471  0.695081  0.98900
 #> 
 #> Resulting formation fractions:
 #>                ff
-#> DMTA_M23  0.11033
-#> DMTA_M27  0.10218
-#> DMTA_M31  0.08784
-#> DMTA_sink 0.69965
+#> DMTA_M23  0.10668
+#> DMTA_M27  0.09481
+#> DMTA_M31  0.08908
+#> DMTA_sink 0.70943
 #> 
 #> Estimated disappearance times:
 #>       DT50   DT90 DT50back DT50_k1 DT50_k2
-#> DMTA 16.59  57.44    17.29   14.71   24.97
-#> M23  53.74 178.51       NA      NA      NA
-#> M27  39.07 129.78       NA      NA      NA
-#> M31  44.06 146.36       NA      NA      NA
plot(f_dmta_nlmixr_saem)
+#> DMTA 12.57  45.43    13.67   11.72   53.73
+#> M23  56.63 188.11       NA      NA      NA
+#> M27  41.29 137.18       NA      NA      NA
+#> M31  43.05 143.01       NA      NA      NA
plot(f_dmta_nlmixr_saem)
 
# }
-- cgit v1.2.1 From ff83d8b2ba623513d92ac90fac4a1b0ec98c2cb5 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 5 Oct 2021 17:33:52 +0200 Subject: Update docs --- docs/dev/reference/dimethenamid_2018.html | 463 +++++++----------------------- 1 file changed, 108 insertions(+), 355 deletions(-) (limited to 'docs/dev/reference/dimethenamid_2018.html') diff --git a/docs/dev/reference/dimethenamid_2018.html b/docs/dev/reference/dimethenamid_2018.html index 919e9363..60c15ade 100644 --- a/docs/dev/reference/dimethenamid_2018.html +++ b/docs/dev/reference/dimethenamid_2018.html @@ -222,7 +222,7 @@ specific pieces of information in the comments.

list("DFOP-SFO3+" = dfop_sfo3_plus), dmta_ds, quiet = TRUE, error_model = "tc") nlmixr_model(f_dmta_mkin_tc) -
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> function () +
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> Warning: number of items to replace is not a multiple of replacement length
#> function () #> { #> ini({ #> DMTA_0 = 98.7132391714013 @@ -263,9 +263,9 @@ specific pieces of information in the comments.

#> k1 = exp(log_k1 + eta.log_k1) #> k2 = exp(log_k2 + eta.log_k2) #> g = expit(g_qlogis + eta.g_qlogis) -#> f_DMTA_tffm0_1 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis) -#> f_DMTA_tffm0_2 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis) -#> f_DMTA_tffm0_3 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis) +#> f_DMTA_to_M23 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis) +#> f_DMTA_to_M23 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis) +#> f_DMTA_to_M23 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis) #> f_DMTA_to_M23 = f_DMTA_tffm0_1 #> f_DMTA_to_M27 = f_DMTA_tffm0_2 * (1 - f_DMTA_tffm0_1) #> f_DMTA_to_M31 = f_DMTA_tffm0_3 * (1 - f_DMTA_tffm0_2) * @@ -289,205 +289,122 @@ specific pieces of information in the comments.

#> M31 ~ add(sigma_low_M31) + prop(rsd_high_M31) #> }) #> } -#> <environment: 0x555559e97ac0>
# The focei fit takes about four minutes on my system +#> <environment: 0x555559d89920>
# 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)) ) -
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 +
#> Warning: number of items to replace is not a multiple of replacement length
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02 #>
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04 #>
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01 -#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:08 +#>
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:09 #>
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 -#>
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:07 +#>
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:06 #>
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #>
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 -#>
#> → compiling inner model...
#>
#> ✔ done
#> → finding duplicate expressions in FD model...
#>
#> → optimizing duplicate expressions in FD model...
#>
#> → compiling EBE model...
#>
#> ✔ done
#> → compiling events FD model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> RxODE 1.1.1 using 8 threads (see ?getRxThreads) -#> no cache: create with `rxCreateCache()`
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation -#> F: Forward difference gradient approximation -#> C: Central difference gradient approximation -#> M: Mixed forward and central difference gradient approximation -#> Unscaled parameters for Omegas=chol(solve(omega)); -#> Diagonals are transformed, as specified by foceiControl(diagXform=) -#> |-----+---------------+-----------+-----------+-----------+-----------| -#> | #| Objective Fun | DMTA_0 | log_k_M23 | log_k_M27 | log_k_M31 | -#> |.....................| log_k1 | log_k2 | g_qlogis |f_DMTA_tffm0_1_qlogis | -#> |.....................|f_DMTA_tffm0_2_qlogis |f_DMTA_tffm0_3_qlogis | sigma_low | rsd_high | -#> |.....................| o1 | o2 | o3 | o4 | -#> |.....................| o5 | o6 | o7 | o8 | -#> |.....................| o9 | o10 |...........|...........| -#> calculating covariance matrix -#> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: last objective function was not at minimum, possible problems in optimization
#> Warning: S matrix non-positive definite
#> Warning: using R matrix to calculate covariance
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
#> user system elapsed -#> 230.015 8.962 238.957
summary(f_dmta_nlmixr_focei) -
#> nlmixr version used for fitting: 2.0.5 -#> mkin version used for pre-fitting: 1.1.0 -#> R version used for fitting: 4.1.1 -#> Date of fit: Thu Sep 16 14:06:55 2021 -#> Date of summary: Thu Sep 16 14:06:55 2021 -#> -#> Equations: -#> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * -#> time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) -#> * DMTA -#> d_M23/dt = + f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M23 * M23 -#> d_M27/dt = + f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M27 * M27 + k_M31 * M31 -#> d_M31/dt = + f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M31 * M31 -#> -#> Data: -#> 563 observations of 4 variable(s) grouped in 6 datasets -#> -#> Degradation model predictions using RxODE -#> -#> Fitted in 238.792 s -#> -#> Variance model: Two-component variance function -#> -#> Mean of starting values for individual parameters: -#> DMTA_0 log_k_M23 log_k_M27 log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 -#> 98.7132 -3.9216 -4.3306 -4.2442 0.1376 0.1388 -#> f_DMTA_ilr_3 log_k1 log_k2 g_qlogis -#> -1.7554 -2.2352 -3.7758 0.4363 -#> -#> Mean of starting values for error model parameters: -#> sigma_low rsd_high -#> 0.70012 0.02577 -#> -#> Fixed degradation parameter values: -#> None -#> -#> Results: -#> -#> Likelihood calculated by focei -#> AIC BIC logLik -#> 1918 2014 -937.2 -#> -#> Optimised parameters: -#> est. lower upper -#> DMTA_0 98.7132 98.6801 98.7464 -#> log_k_M23 -3.9216 -3.9235 -3.9198 -#> log_k_M27 -4.3306 -4.3326 -4.3286 -#> log_k_M31 -4.2442 -4.2461 -4.2422 -#> log_k1 -2.2352 -2.2364 -2.2340 -#> log_k2 -3.7758 -3.7776 -3.7740 -#> g_qlogis 0.4363 0.4358 0.4368 -#> f_DMTA_tffm0_1_qlogis -2.0915 -2.0926 -2.0903 -#> f_DMTA_tffm0_2_qlogis -2.1788 -2.1800 -2.1776 -#> f_DMTA_tffm0_3_qlogis -2.1404 -2.1415 -2.1392 -#> -#> Correlation: -#> DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs -#> log_k_M23 0 -#> log_k_M27 0 0 -#> log_k_M31 0 0 0 -#> log_k1 0 0 0 0 -#> log_k2 0 0 0 0 0 -#> g_qlogis 0 0 0 0 0 0 -#> f_DMTA_tffm0_1_qlogis 0 0 0 0 0 0 0 -#> f_DMTA_tffm0_2_qlogis 0 0 0 0 0 0 0 -#> f_DMTA_tffm0_3_qlogis 0 0 0 0 0 0 0 -#> f_DMTA_0_1 f_DMTA_0_2 -#> log_k_M23 -#> log_k_M27 -#> log_k_M31 -#> log_k1 -#> log_k2 -#> g_qlogis -#> f_DMTA_tffm0_1_qlogis -#> f_DMTA_tffm0_2_qlogis 0 -#> f_DMTA_tffm0_3_qlogis 0 0 -#> -#> Random effects (omega): -#> eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31 -#> eta.DMTA_0 2.327 0.0000 0.0000 0.0000 -#> eta.log_k_M23 0.000 0.5493 0.0000 0.0000 -#> eta.log_k_M27 0.000 0.0000 0.8552 0.0000 -#> eta.log_k_M31 0.000 0.0000 0.0000 0.7457 -#> eta.log_k1 0.000 0.0000 0.0000 0.0000 -#> eta.log_k2 0.000 0.0000 0.0000 0.0000 -#> eta.g_qlogis 0.000 0.0000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.0000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.0000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.0000 0.0000 0.0000 -#> eta.log_k1 eta.log_k2 eta.g_qlogis -#> eta.DMTA_0 0.000 0.000 0.000 -#> eta.log_k_M23 0.000 0.000 0.000 -#> eta.log_k_M27 0.000 0.000 0.000 -#> eta.log_k_M31 0.000 0.000 0.000 -#> eta.log_k1 0.901 0.000 0.000 -#> eta.log_k2 0.000 1.577 0.000 -#> eta.g_qlogis 0.000 0.000 3.102 -#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.000 0.000 -#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.000 0.000 -#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.000 0.000 -#> eta.f_DMTA_tffm0_1_qlogis eta.f_DMTA_tffm0_2_qlogis -#> eta.DMTA_0 0.0 0.0 -#> eta.log_k_M23 0.0 0.0 -#> eta.log_k_M27 0.0 0.0 -#> eta.log_k_M31 0.0 0.0 -#> eta.log_k1 0.0 0.0 -#> eta.log_k2 0.0 0.0 -#> eta.g_qlogis 0.0 0.0 -#> eta.f_DMTA_tffm0_1_qlogis 0.3 0.0 -#> eta.f_DMTA_tffm0_2_qlogis 0.0 0.3 -#> eta.f_DMTA_tffm0_3_qlogis 0.0 0.0 -#> eta.f_DMTA_tffm0_3_qlogis -#> eta.DMTA_0 0.0 -#> eta.log_k_M23 0.0 -#> eta.log_k_M27 0.0 -#> eta.log_k_M31 0.0 -#> eta.log_k1 0.0 -#> eta.log_k2 0.0 -#> eta.g_qlogis 0.0 -#> eta.f_DMTA_tffm0_1_qlogis 0.0 -#> eta.f_DMTA_tffm0_2_qlogis 0.0 -#> eta.f_DMTA_tffm0_3_qlogis 0.3 -#> -#> Variance model: -#> sigma_low rsd_high -#> 0.70012 0.02577 -#> -#> Backtransformed parameters: -#> est. lower upper -#> DMTA_0 98.71324 98.68012 98.74636 -#> k_M23 0.01981 0.01977 0.01985 -#> k_M27 0.01316 0.01313 0.01319 -#> k_M31 0.01435 0.01432 0.01438 -#> f_DMTA_to_M23 0.10993 NA NA -#> f_DMTA_to_M27 0.09049 NA NA -#> f_DMTA_to_M31 0.08414 NA NA -#> k1 0.10698 0.10685 0.10710 -#> k2 0.02292 0.02288 0.02296 -#> g 0.60738 0.60725 0.60751 -#> -#> Resulting formation fractions: -#> ff -#> DMTA_M23 0.10993 -#> DMTA_M27 0.09049 -#> DMTA_M31 0.08414 -#> DMTA_sink 0.71543 -#> -#> Estimated disappearance times: -#> DT50 DT90 DT50back DT50_k1 DT50_k2 -#> DMTA 10.72 60.1 18.09 6.48 30.24 -#> M23 34.99 116.2 NA NA NA -#> M27 52.67 175.0 NA NA NA -#> M31 48.31 160.5 NA NA NA
plot(f_dmta_nlmixr_focei) -
# Using saemix takes about 18 minutes +#>
#> → compiling inner model...
#>
#> ✔ done
#> → finding duplicate expressions in FD model...
#>
#> → optimizing duplicate expressions in FD model...
#>
#> → compiling EBE model...
#>
#> ✔ done
#> → compiling events FD model...
#>
#> ✔ done
#> Model:
#> cmt(DMTA); +#> cmt(M23); +#> cmt(M27); +#> cmt(M31); +#> rx_expr_14~ETA[1]+THETA[1]; +#> DMTA(0)=rx_expr_14; +#> rx_expr_15~ETA[5]+THETA[5]; +#> rx_expr_16~ETA[7]+THETA[7]; +#> rx_expr_17~ETA[6]+THETA[6]; +#> rx_expr_24~exp(rx_expr_15); +#> rx_expr_25~exp(rx_expr_17); +#> rx_expr_29~t*rx_expr_24; +#> rx_expr_30~t*rx_expr_25; +#> rx_expr_31~exp(-(rx_expr_16)); +#> rx_expr_35~1+rx_expr_31; +#> rx_expr_40~1/(rx_expr_35); +#> rx_expr_42~(rx_expr_40); +#> rx_expr_43~1-rx_expr_42; +#> d/dt(DMTA)=-DMTA*(exp(rx_expr_15-rx_expr_29)/(rx_expr_35)+exp(rx_expr_17-rx_expr_30)*(rx_expr_43))/(exp(-t*rx_expr_24)/(rx_expr_35)+exp(-t*rx_expr_25)*(rx_expr_43)); +#> rx_expr_18~ETA[2]+THETA[2]; +#> rx_expr_26~exp(rx_expr_18); +#> d/dt(M23)=-rx_expr_26*M23+DMTA*(exp(rx_expr_15-rx_expr_29)/(rx_expr_35)+exp(rx_expr_17-rx_expr_30)*(rx_expr_43))*f_DMTA_tffm0_1/(exp(-t*rx_expr_24)/(rx_expr_35)+exp(-t*rx_expr_25)*(rx_expr_43)); +#> rx_expr_19~ETA[3]+THETA[3]; +#> rx_expr_20~ETA[4]+THETA[4]; +#> rx_expr_21~1-f_DMTA_tffm0_1; +#> rx_expr_27~exp(rx_expr_19); +#> rx_expr_28~exp(rx_expr_20); +#> d/dt(M27)=-rx_expr_27*M27+rx_expr_28*M31+DMTA*(rx_expr_21)*(exp(rx_expr_15-rx_expr_29)/(rx_expr_35)+exp(rx_expr_17-rx_expr_30)*(rx_expr_43))*f_DMTA_tffm0_2/(exp(-t*rx_expr_24)/(rx_expr_35)+exp(-t*rx_expr_25)*(rx_expr_43)); +#> rx_expr_22~1-f_DMTA_tffm0_2; +#> d/dt(M31)=-rx_expr_28*M31+DMTA*(rx_expr_22)*(rx_expr_21)*(exp(rx_expr_15-rx_expr_29)/(rx_expr_35)+exp(rx_expr_17-rx_expr_30)*(rx_expr_43))*f_DMTA_tffm0_3/(exp(-t*rx_expr_24)/(rx_expr_35)+exp(-t*rx_expr_25)*(rx_expr_43)); +#> rx_expr_0~CMT==4; +#> rx_expr_1~CMT==2; +#> rx_expr_2~CMT==1; +#> rx_expr_3~CMT==3; +#> rx_expr_4~1-(rx_expr_0); +#> rx_expr_5~1-(rx_expr_1); +#> rx_expr_6~1-(rx_expr_3); +#> rx_yj_~(rx_expr_4)*((2*(rx_expr_5)*(rx_expr_2)+2*(rx_expr_1))*(rx_expr_6)+2*(rx_expr_3))+2*(rx_expr_0); +#> rx_expr_7~(rx_expr_1); +#> rx_expr_8~(rx_expr_3); +#> rx_expr_9~(rx_expr_0); +#> rx_expr_13~(rx_expr_5); +#> rx_expr_32~rx_expr_13*(rx_expr_2); +#> rx_lambda_~(rx_expr_4)*((rx_expr_32+rx_expr_7)*(rx_expr_6)+rx_expr_8)+rx_expr_9; +#> rx_hi_~(rx_expr_4)*((rx_expr_32+rx_expr_7)*(rx_expr_6)+rx_expr_8)+rx_expr_9; +#> rx_low_~0; +#> rx_expr_10~M31*(rx_expr_0); +#> rx_expr_11~M27*(rx_expr_3); +#> rx_expr_12~M23*(rx_expr_1); +#> rx_expr_23~DMTA*(rx_expr_5); +#> rx_expr_36~rx_expr_23*(rx_expr_2); +#> rx_pred_=(rx_expr_4)*((rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_3)+((rx_expr_1)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))+(rx_expr_5)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_2))*(rx_expr_6))+(rx_expr_0)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6))); +#> rx_expr_33~Rx_pow_di(THETA[12],2); +#> rx_expr_34~Rx_pow_di(THETA[11],2); +#> rx_r_=(rx_expr_4)*((rx_expr_33*Rx_pow_di(((rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_3)+((rx_expr_1)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))+(rx_expr_5)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_2))*(rx_expr_6)),2)+rx_expr_34)*(rx_expr_3)+((rx_expr_1)*(rx_expr_33*Rx_pow_di(((rx_expr_1)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))+(rx_expr_5)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_2)),2)+rx_expr_34)+(rx_expr_33*Rx_pow_di(((rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_2)),2)+rx_expr_34)*(rx_expr_5)*(rx_expr_2))*(rx_expr_6))+(rx_expr_0)*(rx_expr_33*Rx_pow_di(((rx_expr_4)*((rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_3)+((rx_expr_1)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))+(rx_expr_5)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))*(rx_expr_2))*(rx_expr_6))+(rx_expr_0)*(rx_expr_10+(rx_expr_4)*(rx_expr_11+(rx_expr_12+rx_expr_36)*(rx_expr_6)))),2)+rx_expr_34); +#> DMTA_0=THETA[1]; +#> log_k_M23=THETA[2]; +#> log_k_M27=THETA[3]; +#> log_k_M31=THETA[4]; +#> log_k1=THETA[5]; +#> log_k2=THETA[6]; +#> g_qlogis=THETA[7]; +#> f_DMTA_tffm0_1_qlogis=THETA[8]; +#> f_DMTA_tffm0_2_qlogis=THETA[9]; +#> f_DMTA_tffm0_3_qlogis=THETA[10]; +#> sigma_low=THETA[11]; +#> rsd_high=THETA[12]; +#> eta.DMTA_0=ETA[1]; +#> eta.log_k_M23=ETA[2]; +#> eta.log_k_M27=ETA[3]; +#> eta.log_k_M31=ETA[4]; +#> eta.log_k1=ETA[5]; +#> eta.log_k2=ETA[6]; +#> eta.g_qlogis=ETA[7]; +#> eta.f_DMTA_tffm0_1_qlogis=ETA[8]; +#> eta.f_DMTA_tffm0_2_qlogis=ETA[9]; +#> eta.f_DMTA_tffm0_3_qlogis=ETA[10]; +#> DMTA_0_model=rx_expr_14; +#> k_M23=rx_expr_26; +#> k_M27=rx_expr_27; +#> k_M31=rx_expr_28; +#> k1=rx_expr_24; +#> k2=rx_expr_25; +#> g=1/(rx_expr_35); +#> f_DMTA_to_M23=1/(1+exp(-(ETA[8]+THETA[8]))); +#> f_DMTA_to_M23=1/(1+exp(-(ETA[9]+THETA[9]))); +#> f_DMTA_to_M23=1/(1+exp(-(ETA[10]+THETA[10]))); +#> f_DMTA_to_M23=f_DMTA_tffm0_1; +#> f_DMTA_to_M27=(rx_expr_21)*f_DMTA_tffm0_2; +#> f_DMTA_to_M31=(rx_expr_22)*(rx_expr_21)*f_DMTA_tffm0_3; +#> tad=tad(); +#> dosenum=dosenum();
#> Needed Covariates:
#> [1] "f_DMTA_tffm0_1" "f_DMTA_tffm0_2" "f_DMTA_tffm0_3" "CMT"
#> Error in (function (data, inits, PKpars, model = NULL, pred = NULL, err = NULL, lower = -Inf, upper = Inf, fixed = NULL, skipCov = NULL, control = foceiControl(), thetaNames = NULL, etaNames = NULL, etaMat = NULL, ..., env = NULL, keep = NULL, drop = NULL) { set.seed(control$seed) .pt <- proc.time() RxODE::.setWarnIdSort(FALSE) on.exit(RxODE::.setWarnIdSort(TRUE)) loadNamespace("n1qn1") if (!RxODE::rxIs(control, "foceiControl")) { control <- do.call(foceiControl, control) } if (is.null(env)) { .ret <- new.env(parent = emptyenv()) } else { .ret <- env } .ret$origData <- data .ret$etaNames <- etaNames .ret$thetaFixed <- fixed .ret$control <- control .ret$control$focei.mu.ref <- integer(0) if (is(model, "RxODE") || is(model, "character")) { .ret$ODEmodel <- TRUE if (class(pred) != "function") { stop("pred must be a function specifying the prediction variables in this model.") } } else { .ret$ODEmodel <- TRUE model <- RxODE::rxGetLin(PKpars) pred <- eval(parse(text = "function(){return(Central);}")) } .square <- function(x) x * x .ret$diagXformInv <- c(sqrt = ".square", log = "exp", identity = "identity")[control$diagXform] if (is.null(err)) { err <- eval(parse(text = paste0("function(){err", paste(inits$ERROR[[1]], collapse = ""), "}"))) } .covNames <- .parNames <- c() .ret$adjLik <- control$adjLik .mixed <- !is.null(inits$OMGA) && length(inits$OMGA) > 0 if (!exists("noLik", envir = .ret)) { .atol <- rep(control$atol, length(RxODE::rxModelVars(model)$state)) .rtol <- rep(control$rtol, length(RxODE::rxModelVars(model)$state)) .ssAtol <- rep(control$ssAtol, length(RxODE::rxModelVars(model)$state)) .ssRtol <- rep(control$ssRtol, length(RxODE::rxModelVars(model)$state)) .ret$model <- RxODE::rxSymPySetupPred(model, pred, PKpars, err, grad = (control$derivMethod == 2L), pred.minus.dv = TRUE, sum.prod = control$sumProd, theta.derivs = FALSE, optExpression = control$optExpression, interaction = (control$interaction == 1L), only.numeric = !.mixed, run.internal = TRUE, addProp = control$addProp) if (!is.null(.ret$model$inner)) { .atol <- c(.atol, rep(control$atolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.atol))) .rtol <- c(.rtol, rep(control$rtolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.rtol))) .ret$control$rxControl$atol <- .atol .ret$control$rxControl$rtol <- .rtol .ssAtol <- c(.ssAtol, rep(control$ssAtolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.ssAtol))) .ssRtol <- c(.ssRtol, rep(control$ssRtolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.ssRtol))) .ret$control$rxControl$ssAtol <- .ssAtol .ret$control$rxControl$ssRtol <- .ssRtol } .covNames <- .parNames <- RxODE::rxParams(.ret$model$pred.only) .covNames <- .covNames[regexpr(rex::rex(start, or("THETA", "ETA"), "[", numbers, "]", end), .covNames) == -1] colnames(data) <- sapply(names(data), function(x) { if (any(x == .covNames)) { return(x) } else { return(toupper(x)) } }) .lhs <- c(names(RxODE::rxInits(.ret$model$pred.only)), RxODE::rxLhs(.ret$model$pred.only)) if (length(.lhs) > 0) { .covNames <- .covNames[regexpr(rex::rex(start, or(.lhs), end), .covNames) == -1] } if (length(.covNames) > 0) { if (!all(.covNames %in% names(data))) { message("Model:") RxODE::rxCat(.ret$model$pred.only) message("Needed Covariates:") nlmixrPrint(.covNames) stop("Not all the covariates are in the dataset.") } message("Needed Covariates:") print(.covNames) } .extraPars <- .ret$model$extra.pars } else { if (.ret$noLik) { .atol <- rep(control$atol, length(RxODE::rxModelVars(model)$state)) .rtol <- rep(control$rtol, length(RxODE::rxModelVars(model)$state)) .ret$model <- RxODE::rxSymPySetupPred(model, pred, PKpars, err, grad = FALSE, pred.minus.dv = TRUE, sum.prod = control$sumProd, theta.derivs = FALSE, optExpression = control$optExpression, run.internal = TRUE, only.numeric = TRUE, addProp = control$addProp) if (!is.null(.ret$model$inner)) { .atol <- c(.atol, rep(control$atolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.atol))) .rtol <- c(.rtol, rep(control$rtolSens, length(RxODE::rxModelVars(.ret$model$inner)$state) - length(.rtol))) .ret$control$rxControl$atol <- .atol .ret$control$rxControl$rtol <- .rtol } .covNames <- .parNames <- RxODE::rxParams(.ret$model$pred.only) .covNames <- .covNames[regexpr(rex::rex(start, or("THETA", "ETA"), "[", numbers, "]", end), .covNames) == -1] colnames(data) <- sapply(names(data), function(x) { if (any(x == .covNames)) { return(x) } else { return(toupper(x)) } }) .lhs <- c(names(RxODE::rxInits(.ret$model$pred.only)), RxODE::rxLhs(.ret$model$pred.only)) if (length(.lhs) > 0) { .covNames <- .covNames[regexpr(rex::rex(start, or(.lhs), end), .covNames) == -1] } if (length(.covNames) > 0) { if (!all(.covNames %in% names(data))) { message("Model:") RxODE::rxCat(.ret$model$pred.only) message("Needed Covariates:") nlmixrPrint(.covNames) stop("Not all the covariates are in the dataset.") } message("Needed Covariates:") print(.covNames) } .extraPars <- .ret$model$extra.pars } else { .extraPars <- NULL } } .ret$skipCov <- skipCov if (is.null(skipCov)) { if (is.null(fixed)) { .tmp <- rep(FALSE, length(inits$THTA)) } else { if (length(fixed) < length(inits$THTA)) { .tmp <- c(fixed, rep(FALSE, length(inits$THTA) - length(fixed))) } else { .tmp <- fixed[1:length(inits$THTA)] } } if (exists("uif", envir = .ret)) { .uifErr <- .ret$uif$ini$err[!is.na(.ret$uif$ini$ntheta)] .uifErr <- sapply(.uifErr, function(x) { if (is.na(x)) { return(FALSE) } return(!any(x == c("pow2", "tbs", "tbsYj"))) }) .tmp <- (.tmp | .uifErr) } .ret$skipCov <- c(.tmp, rep(TRUE, length(.extraPars))) .ret$control$focei.mu.ref <- .ret$uif$focei.mu.ref } if (is.null(.extraPars)) { .nms <- c(sprintf("THETA[%s]", seq_along(inits$THTA))) } else { .nms <- c(sprintf("THETA[%s]", seq_along(inits$THTA)), sprintf("ERR[%s]", seq_along(.extraPars))) } if (!is.null(thetaNames) && (length(inits$THTA) + length(.extraPars)) == length(thetaNames)) { .nms <- thetaNames } .ret$thetaNames <- .nms .thetaReset$thetaNames <- .nms if (length(lower) == 1) { lower <- rep(lower, length(inits$THTA)) } else if (length(lower) != length(inits$THTA)) { print(inits$THTA) print(lower) stop("Lower must be a single constant for all the THETA lower bounds, or match the dimension of THETA.") } if (length(upper) == 1) { upper <- rep(upper, length(inits$THTA)) } else if (length(lower) != length(inits$THTA)) { stop("Upper must be a single constant for all the THETA lower bounds, or match the dimension of THETA.") } if (!is.null(.extraPars)) { .ret$model$extra.pars <- eval(call(control$diagXform, .ret$model$extra.pars)) if (length(.ret$model$extra.pars) > 0) { inits$THTA <- c(inits$THTA, .ret$model$extra.pars) .lowerErr <- rep(control$atol[1] * 10, length(.ret$model$extra.pars)) .upperErr <- rep(Inf, length(.ret$model$extra.pars)) lower <- c(lower, .lowerErr) upper <- c(upper, .upperErr) } } if (is.null(data$ID)) stop("\"ID\" not found in data") if (is.null(data$DV)) stop("\"DV\" not found in data") if (is.null(data$EVID)) data$EVID <- 0 if (is.null(data$AMT)) data$AMT <- 0 for (.v in c("TIME", "AMT", "DV", .covNames)) { data[[.v]] <- as.double(data[[.v]]) } .ret$dataSav <- data .ds <- data[data$EVID != 0 & data$EVID != 2, c("ID", "TIME", "AMT", "EVID", .covNames)] .w <- which(tolower(names(data)) == "limit") .limitName <- NULL if (length(.w) == 1L) { .limitName <- names(data)[.w] } .censName <- NULL .w <- which(tolower(names(data)) == "cens") if (length(.w) == 1L) { .censName <- names(data[.w]) } data <- data[data$EVID == 0 | data$EVID == 2, c("ID", "TIME", "DV", "EVID", .covNames, .limitName, .censName)] .w <- which(!(names(.ret$dataSav) %in% c(.covNames, keep))) names(.ret$dataSav)[.w] <- tolower(names(.ret$dataSav[.w])) if (.mixed) { .lh <- .parseOM(inits$OMGA) .nlh <- sapply(.lh, length) .osplt <- rep(1:length(.lh), .nlh) .lini <- list(inits$THTA, unlist(.lh)) .nlini <- sapply(.lini, length) .nsplt <- rep(1:length(.lini), .nlini) .om0 <- .genOM(.lh) if (length(etaNames) == dim(.om0)[1]) { .ret$etaNames <- .ret$etaNames } else { .ret$etaNames <- sprintf("ETA[%d]", seq(1, dim(.om0)[1])) } .ret$rxInv <- RxODE::rxSymInvCholCreate(mat = .om0, diag.xform = control$diagXform) .ret$xType <- .ret$rxInv$xType .om0a <- .om0 .om0a <- .om0a/control$diagOmegaBoundLower .om0b <- .om0 .om0b <- .om0b * control$diagOmegaBoundUpper .om0a <- RxODE::rxSymInvCholCreate(mat = .om0a, diag.xform = control$diagXform) .om0b <- RxODE::rxSymInvCholCreate(mat = .om0b, diag.xform = control$diagXform) .omdf <- data.frame(a = .om0a$theta, m = .ret$rxInv$theta, b = .om0b$theta, diag = .om0a$theta.diag) .omdf$lower <- with(.omdf, ifelse(a > b, b, a)) .omdf$lower <- with(.omdf, ifelse(lower == m, -Inf, lower)) .omdf$lower <- with(.omdf, ifelse(!diag, -Inf, lower)) .omdf$upper <- with(.omdf, ifelse(a < b, b, a)) .omdf$upper <- with(.omdf, ifelse(upper == m, Inf, upper)) .omdf$upper <- with(.omdf, ifelse(!diag, Inf, upper)) .ret$control$nomega <- length(.omdf$lower) .ret$control$neta <- sum(.omdf$diag) .ret$control$ntheta <- length(lower) .ret$control$nfixed <- sum(fixed) lower <- c(lower, .omdf$lower) upper <- c(upper, .omdf$upper) } else { .ret$control$nomega <- 0 .ret$control$neta <- 0 .ret$xType <- -1 .ret$control$ntheta <- length(lower) .ret$control$nfixed <- sum(fixed) } .ret$lower <- lower .ret$upper <- upper .ret$thetaIni <- inits$THTA .scaleC <- double(length(lower)) if (is.null(control$scaleC)) { .scaleC <- rep(NA_real_, length(lower)) } else { .scaleC <- as.double(control$scaleC) if (length(lower) > length(.scaleC)) { .scaleC <- c(.scaleC, rep(NA_real_, length(lower) - length(.scaleC))) } else if (length(lower) < length(.scaleC)) { .scaleC <- .scaleC[seq(1, length(lower))] warning("scaleC control option has more options than estimated population parameters, please check.") } } .ret$scaleC <- .scaleC if (exists("uif", envir = .ret)) { .ini <- as.data.frame(.ret$uif$ini)[!is.na(.ret$uif$ini$err), c("est", "err", "ntheta")] for (.i in seq_along(.ini$err)) { if (is.na(.ret$scaleC[.ini$ntheta[.i]])) { if (any(.ini$err[.i] == c("boxCox", "yeoJohnson", "pow2", "tbs", "tbsYj"))) { .ret$scaleC[.ini$ntheta[.i]] <- 1 } else if (any(.ini$err[.i] == c("prop", "add", "norm", "dnorm", "logn", "dlogn", "lnorm", "dlnorm"))) { .ret$scaleC[.ini$ntheta[.i]] <- 0.5 * abs(.ini$est[.i]) } } } for (.i in .ini$model$extraProps$powTheta) { if (is.na(.ret$scaleC[.i])) .ret$scaleC[.i] <- 1 } .ini <- as.data.frame(.ret$uif$ini) for (.i in .ini$model$extraProps$factorial) { if (is.na(.ret$scaleC[.i])) .ret$scaleC[.i] <- abs(1/digamma(.ini$est[.i] + 1)) } for (.i in .ini$model$extraProps$gamma) { if (is.na(.ret$scaleC[.i])) .ret$scaleC[.i] <- abs(1/digamma(.ini$est[.i])) } for (.i in .ini$model$extraProps$log) { if (is.na(.ret$scaleC[.i])) .ret$scaleC[.i] <- log(abs(.ini$est[.i])) * abs(.ini$est[.i]) } for (.i in .ret$logitThetas) { .b <- .ret$logitThetasLow[.i] .c <- .ret$logitThetasHi[.i] .a <- .ini$est[.i] if (is.na(.ret$scaleC[.i])) { .ret$scaleC[.i] <- 1 * (-.b + .c) * exp(-.a)/((1 + exp(-.a))^2 * (.b + 1 * (-.b + .c)/(1 + exp(-.a)))) } } } names(.ret$thetaIni) <- sprintf("THETA[%d]", seq_along(.ret$thetaIni)) if (is.null(etaMat) & !is.null(control$etaMat)) { .ret$etaMat <- control$etaMat } else { .ret$etaMat <- etaMat } .ret$setupTime <- (proc.time() - .pt)["elapsed"] if (exists("uif", envir = .ret)) { .tmp <- .ret$uif$logThetasList .ret$logThetas <- .tmp[[1]] .ret$logThetasF <- .tmp[[2]] .tmp <- .ret$uif$logitThetasList .ret$logitThetas <- .tmp[[1]] .ret$logitThetasF <- .tmp[[2]] .tmp <- .ret$uif$logitThetasListLow .ret$logitThetasLow <- .tmp[[1]] .ret$logitThetasLowF <- .tmp[[2]] .tmp <- .ret$uif$logitThetasListHi .ret$logitThetasHi <- .tmp[[1]] .ret$logitThetasHiF <- .tmp[[2]] .tmp <- .ret$uif$probitThetasList .ret$probitThetas <- .tmp[[1]] .ret$probitThetasF <- .tmp[[2]] .tmp <- .ret$uif$probitThetasListLow .ret$probitThetasLow <- .tmp[[1]] .ret$probitThetasLowF <- .tmp[[2]] .tmp <- .ret$uif$probitThetasListHi .ret$probitThetasHi <- .tmp[[1]] .ret$probitThetasHiF <- .tmp[[2]] } else { .ret$logThetasF <- integer(0) .ret$logitThetasF <- integer(0) .ret$logitThetasHiF <- numeric(0) .ret$logitThetasLowF <- numeric(0) .ret$logitThetas <- integer(0) .ret$logitThetasHi <- numeric(0) .ret$logitThetasLow <- numeric(0) .ret$probitThetasF <- integer(0) .ret$probitThetasHiF <- numeric(0) .ret$probitThetasLowF <- numeric(0) .ret$probitThetas <- integer(0) .ret$probitThetasHi <- numeric(0) .ret$probitThetasLow <- numeric(0) } if (exists("noLik", envir = .ret)) { if (!.ret$noLik) { .ret$.params <- c(sprintf("THETA[%d]", seq_along(.ret$thetaIni)), sprintf("ETA[%d]", seq(1, dim(.om0)[1]))) .ret$.thetan <- length(.ret$thetaIni) .ret$nobs <- sum(data$EVID == 0) } } .ret$control$printTop <- TRUE .ret$control$nF <- 0 .est0 <- .ret$thetaIni if (!is.null(.ret$model$pred.nolhs)) { .ret$control$predNeq <- length(.ret$model$pred.nolhs$state) } else { .ret$control$predNeq <- 0L } .fitFun <- function(.ret) { this.env <- environment() assign("err", "theta reset", this.env) while (this.env$err == "theta reset") { assign("err", "", this.env) .ret0 <- tryCatch({ foceiFitCpp_(.ret) }, error = function(e) { if (regexpr("theta reset", e$message) != -1) { assign("zeroOuter", FALSE, this.env) assign("zeroGrad", FALSE, this.env) if (regexpr("theta reset0", e$message) != -1) { assign("zeroGrad", TRUE, this.env) } else if (regexpr("theta resetZ", e$message) != -1) { assign("zeroOuter", TRUE, this.env) } assign("err", "theta reset", this.env) } else { assign("err", e$message, this.env) } }) if (this.env$err == "theta reset") { .nm <- names(.ret$thetaIni) .ret$thetaIni <- setNames(.thetaReset$thetaIni + 0, .nm) .ret$rxInv$theta <- .thetaReset$omegaTheta .ret$control$printTop <- FALSE .ret$etaMat <- .thetaReset$etaMat .ret$control$etaMat <- .thetaReset$etaMat .ret$control$maxInnerIterations <- .thetaReset$maxInnerIterations .ret$control$nF <- .thetaReset$nF .ret$control$gillRetC <- .thetaReset$gillRetC .ret$control$gillRet <- .thetaReset$gillRet .ret$control$gillRet <- .thetaReset$gillRet .ret$control$gillDf <- .thetaReset$gillDf .ret$control$gillDf2 <- .thetaReset$gillDf2 .ret$control$gillErr <- .thetaReset$gillErr .ret$control$rEps <- .thetaReset$rEps .ret$control$aEps <- .thetaReset$aEps .ret$control$rEpsC <- .thetaReset$rEpsC .ret$control$aEpsC <- .thetaReset$aEpsC .ret$control$c1 <- .thetaReset$c1 .ret$control$c2 <- .thetaReset$c2 if (this.env$zeroOuter) { message("Posthoc reset") .ret$control$maxOuterIterations <- 0L } else if (this.env$zeroGrad) { message("Theta reset (zero gradient values); Switch to bobyqa") RxODE::rxReq("minqa") .ret$control$outerOptFun <- .bobyqa .ret$control$outerOpt <- -1L } else { message("Theta reset (ETA drift)") } } } if (this.env$err != "") { stop(this.env$err) } else { return(.ret0) } } .ret0 <- try(.fitFun(.ret)) .n <- 1 while (inherits(.ret0, "try-error") && control$maxOuterIterations != 0 && .n <= control$nRetries) { message(sprintf("Restart %s", .n)) .ret$control$nF <- 0 .estNew <- .est0 + 0.2 * .n * abs(.est0) * stats::runif(length(.est0)) - 0.1 * .n .estNew <- sapply(seq_along(.est0), function(.i) { if (.ret$thetaFixed[.i]) { return(.est0[.i]) } else if (.estNew[.i] < lower[.i]) { return(lower + (.Machine$double.eps)^(1/7)) } else if (.estNew[.i] > upper[.i]) { return(upper - (.Machine$double.eps)^(1/7)) } else { return(.estNew[.i]) } }) .ret$thetaIni <- .estNew .ret0 <- try(.fitFun(.ret)) .n <- .n + 1 } if (inherits(.ret0, "try-error")) stop("Could not fit data.") .ret <- .ret0 if (exists("parHistData", .ret)) { .tmp <- .ret$parHistData .tmp <- .tmp[.tmp$type == "Unscaled", names(.tmp) != "type"] .iter <- .tmp$iter .tmp <- .tmp[, names(.tmp) != "iter"] .ret$parHistStacked <- data.frame(stack(.tmp), iter = .iter) names(.ret$parHistStacked) <- c("val", "par", "iter") .ret$parHist <- data.frame(iter = .iter, .tmp) } if (.mixed) { .etas <- .ret$ranef .thetas <- .ret$fixef .pars <- .Call(`_nlmixr_nlmixrParameters`, .thetas, .etas) .ret$shrink <- .Call(`_nlmixr_calcShrinkOnly`, .ret$omega, .pars$eta.lst, length(.etas$ID)) .updateParFixed(.ret) } else { .updateParFixed(.ret) } if (!exists("table", .ret)) { .ret$table <- tableControl() } if (control$calcTables) { .ret <- addTable(.ret, updateObject = "no", keep = keep, drop = drop, table = .ret$table) } .ret})(data = dat, inits = .FoceiInits, PKpars = .pars, model = .mod, pred = function() { return(nlmixr_pred) }, err = uif$error, lower = uif$focei.lower, upper = uif$focei.upper, fixed = uif$focei.fixed, thetaNames = uif$focei.names, etaNames = uif$eta.names, control = control, env = env, keep = .keep, drop = .drop): Not all the covariates are in the dataset.
#> Timing stopped at: 121.4 8.294 129.7
#> Timing stopped at: 121.5 8.294 129.9
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) )
#> Running main SAEM algorithm -#> [1] "Thu Sep 16 14:06:56 2021" +#> [1] "Tue Oct 5 16:58:50 2021" #> .... #> Minimisation finished -#> [1] "Thu Sep 16 14:25:28 2021"
#> user system elapsed -#> 1176.278 0.021 1176.388
+#> [1] "Tue Oct 5 17:17:24 2021"
#> user system elapsed +#> 1181.365 0.031 1181.470
# 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 @@ -498,174 +415,10 @@ specific pieces of information in the comments.

f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)) ) -
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>
#> → generate SAEM model
#> ✔ done
#> 1: 98.3400 -3.5096 -3.3392 -3.7596 -2.2055 -2.7755 1.0281 -2.7872 -2.7223 -2.8341 2.6422 0.7027 0.8124 0.7085 0.8560 1.4980 3.2777 0.3063 0.2850 0.2850 4.1120 0.3716 4.4582 0.3994 4.4820 0.4025 3.7803 0.5780 -#> 500: 97.8212 -4.4030 -4.0872 -4.1289 -2.8278 -4.3505 2.6614 -2.1252 -2.1308 -2.0749 2.9463 1.2933 0.2802 0.3467 0.4814 0.7877 3.0743 0.1508 0.1523 0.3155 0.9557 0.0333 0.4787 0.1073 0.6826 0.0707 0.7849 0.0356
#> Calculating covariance matrix
#>
#> Calculating -2LL by Gaussian quadrature (nnodes=3,nsd=1.6)
#>
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> done
#> user system elapsed -#> 800.784 3.715 149.687
traceplot(f_dmta_nlmixr_saem$nm) +
#> With est = 'saem', a different error model is required for each observed variableChanging the error model to 'obs_tc' (Two-component error for each observed variable)
#> Warning: number of items to replace is not a multiple of replacement length
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> Error in eval(substitute(expr), data, enclos = parent.frame()): Cannot run SAEM since some of the parameters are not mu-referenced (eta.f_DMTA_tffm0_1_qlogis, eta.f_DMTA_tffm0_2_qlogis, eta.f_DMTA_tffm0_3_qlogis)
#> Timing stopped at: 0.849 0.016 0.864
#> Timing stopped at: 1.041 0.016 1.058
traceplot(f_dmta_nlmixr_saem$nm)
#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
summary(f_dmta_nlmixr_saem) -
#> nlmixr version used for fitting: 2.0.5 -#> mkin version used for pre-fitting: 1.1.0 -#> R version used for fitting: 4.1.1 -#> Date of fit: Thu Sep 16 14:29:02 2021 -#> Date of summary: Thu Sep 16 14:29:02 2021 -#> -#> Equations: -#> d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * -#> time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) -#> * DMTA -#> d_M23/dt = + f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M23 * M23 -#> d_M27/dt = + f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M27 * M27 + k_M31 * M31 -#> d_M31/dt = + f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) -#> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * -#> exp(-k2 * time))) * DMTA - k_M31 * M31 -#> -#> Data: -#> 563 observations of 4 variable(s) grouped in 6 datasets -#> -#> Degradation model predictions using RxODE -#> -#> Fitted in 149.421 s -#> -#> Variance model: Two-component variance function -#> -#> Mean of starting values for individual parameters: -#> DMTA_0 log_k_M23 log_k_M27 log_k_M31 f_DMTA_ilr_1 f_DMTA_ilr_2 -#> 98.7132 -3.9216 -4.3306 -4.2442 0.1376 0.1388 -#> f_DMTA_ilr_3 log_k1 log_k2 g_qlogis -#> -1.7554 -2.2352 -3.7758 0.4363 -#> -#> Mean of starting values for error model parameters: -#> sigma_low_DMTA rsd_high_DMTA sigma_low_M23 rsd_high_M23 sigma_low_M27 -#> 0.70012 0.02577 0.70012 0.02577 0.70012 -#> rsd_high_M27 sigma_low_M31 rsd_high_M31 -#> 0.02577 0.70012 0.02577 -#> -#> Fixed degradation parameter values: -#> None -#> -#> Results: -#> -#> Likelihood calculated by focei -#> AIC BIC logLik -#> 1953 2074 -948.3 -#> -#> Optimised parameters: -#> est. lower upper -#> DMTA_0 97.821 95.862 99.780 -#> log_k_M23 -4.403 -5.376 -3.430 -#> log_k_M27 -4.087 -4.545 -3.629 -#> log_k_M31 -4.129 -4.639 -3.618 -#> log_k1 -2.828 -3.389 -2.266 -#> log_k2 -4.351 -5.472 -3.229 -#> g_qlogis 2.661 0.824 4.499 -#> f_DMTA_tffm0_1_qlogis -2.125 -2.449 -1.801 -#> f_DMTA_tffm0_2_qlogis -2.131 -2.468 -1.794 -#> f_DMTA_tffm0_3_qlogis -2.075 -2.540 -1.610 -#> -#> Correlation: -#> DMTA_0 l__M23 l__M27 l__M31 log_k1 log_k2 g_qlgs -#> log_k_M23 -0.019 -#> log_k_M27 -0.028 0.004 -#> log_k_M31 -0.019 0.003 0.075 -#> log_k1 0.038 -0.004 -0.006 -0.003 -#> log_k2 0.046 0.011 0.008 0.009 0.068 -#> g_qlogis -0.067 0.004 0.006 0.001 -0.076 -0.409 -#> f_DMTA_tffm0_1_qlogis -0.062 0.055 0.006 0.004 -0.008 -0.004 0.012 -#> f_DMTA_tffm0_2_qlogis -0.062 0.010 0.058 -0.034 -0.008 -0.007 0.014 -#> f_DMTA_tffm0_3_qlogis -0.052 0.009 0.056 0.071 -0.006 -0.001 0.008 -#> f_DMTA_0_1 f_DMTA_0_2 -#> log_k_M23 -#> log_k_M27 -#> log_k_M31 -#> log_k1 -#> log_k2 -#> g_qlogis -#> f_DMTA_tffm0_1_qlogis -#> f_DMTA_tffm0_2_qlogis 0.017 -#> f_DMTA_tffm0_3_qlogis 0.014 -0.005 -#> -#> Random effects (omega): -#> eta.DMTA_0 eta.log_k_M23 eta.log_k_M27 eta.log_k_M31 -#> eta.DMTA_0 2.946 0.000 0.0000 0.0000 -#> eta.log_k_M23 0.000 1.293 0.0000 0.0000 -#> eta.log_k_M27 0.000 0.000 0.2802 0.0000 -#> eta.log_k_M31 0.000 0.000 0.0000 0.3467 -#> eta.log_k1 0.000 0.000 0.0000 0.0000 -#> eta.log_k2 0.000 0.000 0.0000 0.0000 -#> eta.g_qlogis 0.000 0.000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_1_qlogis 0.000 0.000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_2_qlogis 0.000 0.000 0.0000 0.0000 -#> eta.f_DMTA_tffm0_3_qlogis 0.000 0.000 0.0000 0.0000 -#> eta.log_k1 eta.log_k2 eta.g_qlogis -#> eta.DMTA_0 0.0000 0.0000 0.000 -#> eta.log_k_M23 0.0000 0.0000 0.000 -#> eta.log_k_M27 0.0000 0.0000 0.000 -#> eta.log_k_M31 0.0000 0.0000 0.000 -#> eta.log_k1 0.4814 0.0000 0.000 -#> eta.log_k2 0.0000 0.7877 0.000 -#> eta.g_qlogis 0.0000 0.0000 3.074 -#> eta.f_DMTA_tffm0_1_qlogis 0.0000 0.0000 0.000 -#> eta.f_DMTA_tffm0_2_qlogis 0.0000 0.0000 0.000 -#> eta.f_DMTA_tffm0_3_qlogis 0.0000 0.0000 0.000 -#> eta.f_DMTA_tffm0_1_qlogis eta.f_DMTA_tffm0_2_qlogis -#> eta.DMTA_0 0.0000 0.0000 -#> eta.log_k_M23 0.0000 0.0000 -#> eta.log_k_M27 0.0000 0.0000 -#> eta.log_k_M31 0.0000 0.0000 -#> eta.log_k1 0.0000 0.0000 -#> eta.log_k2 0.0000 0.0000 -#> eta.g_qlogis 0.0000 0.0000 -#> eta.f_DMTA_tffm0_1_qlogis 0.1508 0.0000 -#> eta.f_DMTA_tffm0_2_qlogis 0.0000 0.1523 -#> eta.f_DMTA_tffm0_3_qlogis 0.0000 0.0000 -#> eta.f_DMTA_tffm0_3_qlogis -#> eta.DMTA_0 0.0000 -#> eta.log_k_M23 0.0000 -#> eta.log_k_M27 0.0000 -#> eta.log_k_M31 0.0000 -#> eta.log_k1 0.0000 -#> eta.log_k2 0.0000 -#> eta.g_qlogis 0.0000 -#> eta.f_DMTA_tffm0_1_qlogis 0.0000 -#> eta.f_DMTA_tffm0_2_qlogis 0.0000 -#> eta.f_DMTA_tffm0_3_qlogis 0.3155 -#> -#> Variance model: -#> sigma_low_DMTA rsd_high_DMTA sigma_low_M23 rsd_high_M23 sigma_low_M27 -#> 0.95572 0.03325 0.47871 0.10733 0.68264 -#> rsd_high_M27 sigma_low_M31 rsd_high_M31 -#> 0.07072 0.78486 0.03557 -#> -#> Backtransformed parameters: -#> est. lower upper -#> DMTA_0 97.82122 95.862233 99.78020 -#> k_M23 0.01224 0.004625 0.03239 -#> k_M27 0.01679 0.010615 0.02654 -#> k_M31 0.01610 0.009664 0.02683 -#> f_DMTA_to_M23 0.10668 NA NA -#> f_DMTA_to_M27 0.09481 NA NA -#> f_DMTA_to_M31 0.08908 NA NA -#> k1 0.05914 0.033731 0.10370 -#> k2 0.01290 0.004204 0.03958 -#> g 0.93471 0.695081 0.98900 -#> -#> Resulting formation fractions: -#> ff -#> DMTA_M23 0.10668 -#> DMTA_M27 0.09481 -#> DMTA_M31 0.08908 -#> DMTA_sink 0.70943 -#> -#> Estimated disappearance times: -#> DT50 DT90 DT50back DT50_k1 DT50_k2 -#> DMTA 12.57 45.43 13.67 11.72 53.73 -#> M23 56.63 188.11 NA NA NA -#> M27 41.29 137.18 NA NA NA -#> M31 43.05 143.01 NA NA NA
plot(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
# }