Aerobic soil degradation data on dimethenamid and dimethenamid-P from the EU assessment in 2018

@@ -157,294 +93,368 @@ context of pesticide registrations, as the use of the data may be constrained by data protection regulations.

dimethenamid_2018

- - -

Format

- -

An mkindsg object grouping seven datasets with some meta information

Source

dimethenamid_2018

Format

An mkindsg object grouping seven datasets with some meta information

Source

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 -https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716

Details

- +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 in the 'dataset_generation' directory. In the code, page numbers are given for specific pieces of information in the comments.

Examples

print(dimethenamid_2018)
-
#> <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 
-#>     3     7     7     7     4 
-#> Time normalisation factors $f_time_norm:
-#> [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 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[["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)
#> Warning: number of items to replace is not a multiple of replacement length
#> function () 
-#> {
-#>     ini({
-#>         DMTA_0 = 98.7132391714013
-#>         eta.DMTA_0 ~ 2.32692496033921
-#>         log_k_M23 = -3.92162409637283
-#>         eta.log_k_M23 ~ 0.549278519419884
-#>         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.17879574608926
-#>         eta.f_DMTA_tffm0_2_qlogis ~ 0.3
-#>         f_DMTA_tffm0_3_qlogis = -2.14036526460782
-#>         eta.f_DMTA_tffm0_3_qlogis ~ 0.3
-#>         sigma_low_DMTA = 0.700117227383809
-#>         rsd_high_DMTA = 0.0257724286053519
-#>         sigma_low_M23 = 0.700117227383809
-#>         rsd_high_M23 = 0.0257724286053519
-#>         sigma_low_M27 = 0.700117227383809
-#>         rsd_high_M27 = 0.0257724286053519
-#>         sigma_low_M31 = 0.700117227383809
-#>         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_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) * 
-#>             (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: 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))
-)
-
#> 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:09 
-#> 
#> → finding 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
#> 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] "Tue Oct  5 16:58:50 2021"
-#> ....
-#>     Minimisation finished
-#> [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
-# 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)
#> 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)
-
#> 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
# }
-

Examples

print(dimethenamid_2018)
+#> <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 
+#>     3     7     7     7     4 
+#> Time normalisation factors $f_time_norm:
+#> [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 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[["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)
+#> Warning: number of items to replace is not a multiple of replacement length
+#> function () 
+#> {
+#>     ini({
+#>         DMTA_0 = 99
+#>         eta.DMTA_0 ~ 2.3
+#>         log_k_M23 = -3.9
+#>         eta.log_k_M23 ~ 0.55
+#>         log_k_M27 = -4.3
+#>         eta.log_k_M27 ~ 0.86
+#>         log_k_M31 = -4.2
+#>         eta.log_k_M31 ~ 0.75
+#>         log_k1 = -2.2
+#>         eta.log_k1 ~ 0.9
+#>         log_k2 = -3.8
+#>         eta.log_k2 ~ 1.6
+#>         g_qlogis = 0.44
+#>         eta.g_qlogis ~ 3.1
+#>         f_DMTA_tffm0_1_qlogis = -2.1
+#>         eta.f_DMTA_tffm0_1_qlogis ~ 0.3
+#>         f_DMTA_tffm0_2_qlogis = -2.2
+#>         eta.f_DMTA_tffm0_2_qlogis ~ 0.3
+#>         f_DMTA_tffm0_3_qlogis = -2.1
+#>         eta.f_DMTA_tffm0_3_qlogis ~ 0.3
+#>         sigma_low_DMTA = 0.7
+#>         rsd_high_DMTA = 0.026
+#>         sigma_low_M23 = 0.7
+#>         rsd_high_M23 = 0.026
+#>         sigma_low_M27 = 0.7
+#>         rsd_high_M27 = 0.026
+#>         sigma_low_M31 = 0.7
+#>         rsd_high_M31 = 0.026
+#>     })
+#>     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_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) * 
+#>             (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: 0x55555fca3790>
+# 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))
+)
+#> 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:01 
+#> 
+#> → calculate sensitivities
+#> [====|====|====|====|====|====|====|====|====|====] 0:00:03 
+#> 
+#> → 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: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
+#> 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: 119.8 9.331 129.2
+#> Timing stopped at: 120 9.331 129.3
+summary(f_dmta_nlmixr_focei)
+#> Error in summary(f_dmta_nlmixr_focei): object 'f_dmta_nlmixr_focei' not found
+plot(f_dmta_nlmixr_focei)
+#> Error in plot(f_dmta_nlmixr_focei): object 'f_dmta_nlmixr_focei' not found
+# Using saemix takes about 18 minutes
+system.time(
+  f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE)
+)
+#> DINTDY-  T (=R1) illegal      
+#> In above message, R1 = 115.507
+#>  
+#>       T not in interval TCUR - HU (= R1) to TCUR (=R2)      
+#> In above message, R1 = 112.133, R2 = 113.577
+#>  
+#> DLSODA-  At T (=R1), too much accuracy requested  
+#>       for precision of machine..  See TOLSF (=R2) 
+#> In above message, R1 = 55.3899, R2 = nan
+#>  
+#> Error in out[available, var]: (subscript) logical subscript too long
+#> Timing stopped at: 11.84 0.008 11.85
+#> Timing stopped at: 12.16 0.008 12.17
+
+# 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
+#> 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.892 0.004 0.896
+#> Timing stopped at: 1.096 0.005 1.1
+traceplot(f_dmta_nlmixr_saem$nm)
+#> Error in traceplot(f_dmta_nlmixr_saem$nm): could not find function "traceplot"
+summary(f_dmta_nlmixr_saem)
+#> Error in summary(f_dmta_nlmixr_saem): object 'f_dmta_nlmixr_saem' not found
+plot(f_dmta_nlmixr_saem)
+#> Error in plot(f_dmta_nlmixr_saem): object 'f_dmta_nlmixr_saem' not found
+# }
+