From b5b446b718b15ccaae5b197e147fc1358f0f564e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 6 Nov 2020 00:03:29 +0100 Subject: Fast analytical solutions for saemix, update.mmkin Also, use logit transformation for g and for solitary formation fractions, addressing #10. --- docs/dev/reference/saemix.html | 974 +++++++++++++++++++++++++++++++---------- 1 file changed, 748 insertions(+), 226 deletions(-) (limited to 'docs/dev/reference/saemix.html') diff --git a/docs/dev/reference/saemix.html b/docs/dev/reference/saemix.html index c8cf9fab..5dacefc9 100644 --- a/docs/dev/reference/saemix.html +++ b/docs/dev/reference/saemix.html @@ -190,16 +190,15 @@ mmkin. Starting variances of the random effects (argument omega.init) are the variances of the deviations of the parameters from these mean values.

Examples

- 
ds <- lapply(experimental_data_for_UBA_2019[6:10],
- function(x) subset(x$data[c("name", "time", "value")]))
-names(ds) <- paste("Dataset", 6:10)
-sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
-  A1 = mkinsub("SFO"))
-
#> Successfully compiled differential equation model from auto-generated C code.
# \dontrun{
-f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE)
+    
# \dontrun{
 library(saemix)
 
#> Package saemix, version 3.1.9000
-#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
m_saemix <- saemix_model(f_mmkin, cores = 1)
+#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
ds <- lapply(experimental_data_for_UBA_2019[6:10],
+ function(x) subset(x$data[c("name", "time", "value")]))
+names(ds) <- paste("Dataset", 6:10)
+f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
+  state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
+m_saemix_p0_fixed <- saemix_model(f_mmkin_parent_p0_fixed["FOMC", ])
 
#> 
 #> 
 #> The following SaemixModel object was successfully created:
@@ -208,59 +207,236 @@ variances of the deviations of the parameters from these mean values.

 #>   Model function:  Mixed model generated from mmkin object  Model type:  structural
 #> function (psi, id, xidep) 
 #> {
-#>     uid <- unique(id)
-#>     res_list <- parallel::mclapply(uid, function(i) {
-#>         transparms_optim <- psi[i, ]
-#>         names(transparms_optim) <- names(degparms_optim)
-#>         odeini_optim <- transparms_optim[odeini_optim_parm_names]
-#>         names(odeini_optim) <- gsub("_0$", "", odeini_optim_parm_names)
-#>         odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
-#>         ode_transparms_optim_names <- setdiff(names(transparms_optim), 
-#>             odeini_optim_parm_names)
-#>         odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], 
-#>             mkin_model, transform_rates = object[[1]]$transform_rates, 
-#>             transform_fractions = object[[1]]$transform_fractions)
-#>         odeparms <- c(odeparms_optim, odeparms_fixed)
-#>         xidep_i <- subset(xidep, id == i)
-#>         if (analytical) {
-#>             out_values <- mkin_model$deg_func(xidep_i, odeini, 
-#>                 odeparms)
-#>         }
-#>         else {
-#>             i_time <- xidep_i$time
-#>             i_name <- xidep_i$name
-#>             out_wide <- mkinpredict(mkin_model, odeparms = odeparms, 
-#>                 odeini = odeini, solution_type = object[[1]]$solution_type, 
-#>                 outtimes = sort(unique(i_time)))
-#>             out_index <- cbind(as.character(i_time), as.character(i_name))
-#>             out_values <- out_wide[out_index]
-#>         }
-#>         return(out_values)
-#>     }, mc.cores = cores)
-#>     res <- unlist(res_list)
-#>     return(res)
+#>     odeini_fixed/(xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 
+#>         1])
+#> }
+#> <bytecode: 0x5555599945b8>
+#> <environment: 0x555559984388>
+#>   Nb of parameters: 2 
+#>       parameter names:  log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] log_alpha normal       Estimated
+#> [2,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           log_alpha log_beta
+#> log_alpha         1        0
+#> log_beta          0        1
+#>   Error model: constant , initial values: a.1=2.95893806804889 
+#>     No covariate in the model.
+#>     Initial values
+#>              log_alpha log_beta
+#> Pop.CondInit -0.347996  1.66788
d_saemix_parent <- saemix_data(f_mmkin_parent_p0_fixed)
+
#> 
+#> 
+#> The following SaemixData object was successfully created:
+#> 
+#> Object of class SaemixData
+#>     longitudinal data for use with the SAEM algorithm
+#> Dataset ds_saemix 
+#>     Structured data: value ~ time + name | ds 
+#>     X variable for graphs: time () 
saemix_options <- list(seed = 123456, displayProgress = FALSE,
+  save = FALSE, save.graphs = FALSE, nbiter.saemix = c(200, 80))
+f_saemix_p0_fixed <- saemix(m_saemix_p0_fixed, d_saemix_parent, saemix_options)
+
#> Running main SAEM algorithm
+#> [1] "Thu Nov  5 23:53:29 2020"
+#> ..
+#>     Minimisation finished
+#> [1] "Thu Nov  5 23:53:30 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> -----------------------------------
+#> ----          Data             ----
+#> -----------------------------------
+#> Object of class SaemixData
+#>     longitudinal data for use with the SAEM algorithm
+#> Dataset ds_saemix 
+#>     Structured data: value ~ time + name | ds 
+#>     X variable for graphs: time () 
+#> Dataset characteristics:
+#>     number of subjects:     5 
+#>     number of observations: 90 
+#>     average/min/max nb obs: 18.00  /  16  /  20 
+#> First 10 lines of data:
+#>           ds time   name value mdv cens occ ytype
+#> 1  Dataset 6    0 parent  97.2   0    0   1     1
+#> 2  Dataset 6    0 parent  96.4   0    0   1     1
+#> 3  Dataset 6    3 parent  71.1   0    0   1     1
+#> 4  Dataset 6    3 parent  69.2   0    0   1     1
+#> 5  Dataset 6    6 parent  58.1   0    0   1     1
+#> 6  Dataset 6    6 parent  56.6   0    0   1     1
+#> 7  Dataset 6   10 parent  44.4   0    0   1     1
+#> 8  Dataset 6   10 parent  43.4   0    0   1     1
+#> 9  Dataset 6   20 parent  33.3   0    0   1     1
+#> 10 Dataset 6   20 parent  29.2   0    0   1     1
+#> -----------------------------------
+#> ----          Model            ----
+#> -----------------------------------
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     odeini_fixed/(xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 
+#>         1])
+#> }
+#> <bytecode: 0x5555599945b8>
+#> <environment: 0x555559984388>
+#>   Nb of parameters: 2 
+#>       parameter names:  log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] log_alpha normal       Estimated
+#> [2,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           log_alpha log_beta
+#> log_alpha         1        0
+#> log_beta          0        1
+#>   Error model: constant , initial values: a.1=2.95893806804889 
+#>     No covariate in the model.
+#>     Initial values
+#>              log_alpha log_beta
+#> Pop.CondInit -0.347996  1.66788
+#> -----------------------------------
+#> ----    Key algorithm options  ----
+#> -----------------------------------
+#>     Estimation of individual parameters (MAP)
+#>     Estimation of standard errors and linearised log-likelihood
+#>     Estimation of log-likelihood by importance sampling
+#>     Number of iterations:  K1=200, K2=80 
+#>     Number of chains:  10 
+#>     Seed:  123456 
+#>     Number of MCMC iterations for IS:  5000 
+#>     Simulations:
+#>         nb of simulated datasets used for npde:  1000 
+#>         nb of simulated datasets used for VPC:  100 
+#>     Input/output
+#>         save the results to a file:  FALSE 
+#>         save the graphs to files:  FALSE 
+#> ----------------------------------------------------
+#> ----                  Results                   ----
+#> ----------------------------------------------------
+#> -----------------  Fixed effects  ------------------
+#> ----------------------------------------------------
+#>   Parameter Estimate SE   CV(%)
+#>   log_alpha -0.33    0.30 91.6 
+#>   log_beta   1.70    0.21 12.4 
+#> a a.1        3.15    0.25  7.9 
+#> ----------------------------------------------------
+#> -----------  Variance of random effects  -----------
+#> ----------------------------------------------------
+#>           Parameter        Estimate SE   CV(%)
+#> log_alpha omega2.log_alpha 0.44     0.28 65   
+#> log_beta  omega2.log_beta  0.18     0.14 79   
+#> ----------------------------------------------------
+#> ------  Correlation matrix of random effects  ------
+#> ----------------------------------------------------
+#>                  omega2.log_alpha omega2.log_beta
+#> omega2.log_alpha 1                0              
+#> omega2.log_beta  0                1              
+#> ----------------------------------------------------
+#> ---------------  Statistical criteria  -------------
+#> ----------------------------------------------------
+#> Likelihood computed by linearisation
+#>       -2LL= 501.6082 
+#>       AIC = 511.6082 
+#>       BIC = 509.6554 
+#> 
+#> Likelihood computed by importance sampling
+#>       -2LL= 501.7 
+#>       AIC = 511.7 
+#>       BIC = 509.7472 
+#> ----------------------------------------------------

+f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
+m_saemix_sfo <- saemix_model(f_mmkin_parent["SFO", ])
+
#> 
+#> 
+#> The following SaemixModel object was successfully created:
+#> 
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     psi[id, 1] * exp(-exp(psi[id, 2]) * xidep[, "time"])
 #> }
-#> <bytecode: 0x55555d62aeb8>
-#> <environment: 0x55555e35c170>
+#> <bytecode: 0x55555998d588>
+#> <environment: 0x55555c0f4ae8>
+#>   Nb of parameters: 2 
+#>       parameter names:  parent_0 log_k_parent 
+#>       distribution:
+#>      Parameter    Distribution Estimated
+#> [1,] parent_0     normal       Estimated
+#> [2,] log_k_parent normal       Estimated
+#>   Variance-covariance matrix:
+#>              parent_0 log_k_parent
+#> parent_0            1            0
+#> log_k_parent        0            1
+#>   Error model: constant , initial values: a.1=5.76827561471585 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0 log_k_parent
+#> Pop.CondInit 86.39406    -3.215063
m_saemix_fomc <- saemix_model(f_mmkin_parent["FOMC", ])
+
#> 
+#> 
+#> The following SaemixModel object was successfully created:
+#> 
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     psi[id, 1]/(xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 
+#>         2])
+#> }
+#> <bytecode: 0x55555998dc50>
+#> <environment: 0x5555595d7668>
+#>   Nb of parameters: 3 
+#>       parameter names:  parent_0 log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_alpha normal       Estimated
+#> [3,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           parent_0 log_alpha log_beta
+#> parent_0         1         0        0
+#> log_alpha        0         1        0
+#> log_beta         0         0        1
+#>   Error model: constant , initial values: a.1=1.91976382242696 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0  log_alpha log_beta
+#> Pop.CondInit 94.43884 -0.2226095 2.048192
m_saemix_dfop <- saemix_model(f_mmkin_parent["DFOP", ])
+
#> 
+#> 
+#> The following SaemixModel object was successfully created:
+#> 
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     g <- plogis(psi[id, 4])
+#>     t = xidep[, "time"]
+#>     psi[id, 1] * (g * exp(-exp(psi[id, 2]) * t) + (1 - g) * exp(-exp(psi[id, 
+#>         3]) * t))
+#> }
+#> <bytecode: 0x55555998e548>
+#> <environment: 0x555558225bf0>
 #>   Nb of parameters: 4 
-#>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
+#>       parameter names:  parent_0 log_k1 log_k2 g_qlogis 
 #>       distribution:
-#>      Parameter      Distribution Estimated
-#> [1,] parent_0       normal       Estimated
-#> [2,] log_k_parent   normal       Estimated
-#> [3,] log_k_A1       normal       Estimated
-#> [4,] f_parent_ilr_1 normal       Estimated
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_k1    normal       Estimated
+#> [3,] log_k2    normal       Estimated
+#> [4,] g_qlogis  normal       Estimated
 #>   Variance-covariance matrix:
-#>                parent_0 log_k_parent log_k_A1 f_parent_ilr_1
-#> parent_0              1            0        0              0
-#> log_k_parent          0            1        0              0
-#> log_k_A1              0            0        1              0
-#> f_parent_ilr_1        0            0        0              1
-#>   Error model: constant , initial values: a.1=4.97259024646577 
+#>          parent_0 log_k1 log_k2 g_qlogis
+#> parent_0        1      0      0        0
+#> log_k1          0      1      0        0
+#> log_k2          0      0      1        0
+#> g_qlogis        0      0      0        1
+#>   Error model: constant , initial values: a.1=1.94671278396371 
 #>     No covariate in the model.
 #>     Initial values
-#>              parent_0 log_k_parent  log_k_A1 f_parent_ilr_1
-#> Pop.CondInit 86.53449    -3.207005 -3.060308      -1.920449
d_saemix <- saemix_data(f_mmkin)
+#>              parent_0    log_k1    log_k2 g_qlogis
+#> Pop.CondInit 94.08322 -1.834163 -4.210797  0.11002
d_saemix_parent <- saemix_data(f_mmkin_parent["SFO", ])
 
#> 
 #> 
 #> The following SaemixData object was successfully created:
@@ -269,15 +445,12 @@ variances of the deviations of the parameters from these mean values.

 #>     longitudinal data for use with the SAEM algorithm
 #> Dataset ds_saemix 
 #>     Structured data: value ~ time + name | ds 
-#>     X variable for graphs: time () 
saemix_options <- list(seed = 123456,
-  save = FALSE, save.graphs = FALSE, displayProgress = FALSE,
-  nbiter.saemix = c(200, 80))
-f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
+#>     X variable for graphs: time () 
f_saemix_sfo <- saemix(m_saemix_sfo, d_saemix_parent, saemix_options)
 
#> Running main SAEM algorithm
-#> [1] "Thu Nov  5 08:26:39 2020"
+#> [1] "Thu Nov  5 23:53:31 2020"
 #> ..
 #>     Minimisation finished
-#> [1] "Thu Nov  5 08:28:33 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> [1] "Thu Nov  5 23:53:32 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
 #> -----------------------------------
 #> ----          Data             ----
 #> -----------------------------------
@@ -288,8 +461,8 @@ variances of the deviations of the parameters from these mean values.

 #>     X variable for graphs: time () 
 #> Dataset characteristics:
 #>     number of subjects:     5 
-#>     number of observations: 170 
-#>     average/min/max nb obs: 34.00  /  30  /  38 
+#>     number of observations: 90 
+#>     average/min/max nb obs: 18.00  /  16  /  20 
 #> First 10 lines of data:
 #>           ds time   name value mdv cens occ ytype
 #> 1  Dataset 6    0 parent  97.2   0    0   1     1
@@ -309,59 +482,248 @@ variances of the deviations of the parameters from these mean values.

 #>   Model function:  Mixed model generated from mmkin object  Model type:  structural
 #> function (psi, id, xidep) 
 #> {
-#>     uid <- unique(id)
-#>     res_list <- parallel::mclapply(uid, function(i) {
-#>         transparms_optim <- psi[i, ]
-#>         names(transparms_optim) <- names(degparms_optim)
-#>         odeini_optim <- transparms_optim[odeini_optim_parm_names]
-#>         names(odeini_optim) <- gsub("_0$", "", odeini_optim_parm_names)
-#>         odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
-#>         ode_transparms_optim_names <- setdiff(names(transparms_optim), 
-#>             odeini_optim_parm_names)
-#>         odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], 
-#>             mkin_model, transform_rates = object[[1]]$transform_rates, 
-#>             transform_fractions = object[[1]]$transform_fractions)
-#>         odeparms <- c(odeparms_optim, odeparms_fixed)
-#>         xidep_i <- subset(xidep, id == i)
-#>         if (analytical) {
-#>             out_values <- mkin_model$deg_func(xidep_i, odeini, 
-#>                 odeparms)
-#>         }
-#>         else {
-#>             i_time <- xidep_i$time
-#>             i_name <- xidep_i$name
-#>             out_wide <- mkinpredict(mkin_model, odeparms = odeparms, 
-#>                 odeini = odeini, solution_type = object[[1]]$solution_type, 
-#>                 outtimes = sort(unique(i_time)))
-#>             out_index <- cbind(as.character(i_time), as.character(i_name))
-#>             out_values <- out_wide[out_index]
-#>         }
-#>         return(out_values)
-#>     }, mc.cores = cores)
-#>     res <- unlist(res_list)
-#>     return(res)
+#>     psi[id, 1] * exp(-exp(psi[id, 2]) * xidep[, "time"])
 #> }
-#> <bytecode: 0x55555d62aeb8>
-#> <environment: 0x55555e35c170>
+#> <bytecode: 0x55555998d588>
+#> <environment: 0x55555c0f4ae8>
+#>   Nb of parameters: 2 
+#>       parameter names:  parent_0 log_k_parent 
+#>       distribution:
+#>      Parameter    Distribution Estimated
+#> [1,] parent_0     normal       Estimated
+#> [2,] log_k_parent normal       Estimated
+#>   Variance-covariance matrix:
+#>              parent_0 log_k_parent
+#> parent_0            1            0
+#> log_k_parent        0            1
+#>   Error model: constant , initial values: a.1=5.76827561471585 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0 log_k_parent
+#> Pop.CondInit 86.39406    -3.215063
+#> -----------------------------------
+#> ----    Key algorithm options  ----
+#> -----------------------------------
+#>     Estimation of individual parameters (MAP)
+#>     Estimation of standard errors and linearised log-likelihood
+#>     Estimation of log-likelihood by importance sampling
+#>     Number of iterations:  K1=200, K2=80 
+#>     Number of chains:  10 
+#>     Seed:  123456 
+#>     Number of MCMC iterations for IS:  5000 
+#>     Simulations:
+#>         nb of simulated datasets used for npde:  1000 
+#>         nb of simulated datasets used for VPC:  100 
+#>     Input/output
+#>         save the results to a file:  FALSE 
+#>         save the graphs to files:  FALSE 
+#> ----------------------------------------------------
+#> ----                  Results                   ----
+#> ----------------------------------------------------
+#> -----------------  Fixed effects  ------------------
+#> ----------------------------------------------------
+#>   Parameter    Estimate SE   CV(%)
+#>   parent_0     85.8     1.85  2.2 
+#>   log_k_parent -3.2     0.59 18.3 
+#> a a.1           6.2     0.49  7.9 
+#> ----------------------------------------------------
+#> -----------  Variance of random effects  -----------
+#> ----------------------------------------------------
+#>              Parameter           Estimate SE   CV(%)
+#> parent_0     omega2.parent_0     7.8      10.7 138  
+#> log_k_parent omega2.log_k_parent 1.7       1.1  64  
+#> ----------------------------------------------------
+#> ------  Correlation matrix of random effects  ------
+#> ----------------------------------------------------
+#>                     omega2.parent_0 omega2.log_k_parent
+#> omega2.parent_0     1               0                  
+#> omega2.log_k_parent 0               1                  
+#> ----------------------------------------------------
+#> ---------------  Statistical criteria  -------------
+#> ----------------------------------------------------
+#> Likelihood computed by linearisation
+#>       -2LL= 615.4074 
+#>       AIC = 625.4074 
+#>       BIC = 623.4546 
+#> 
+#> Likelihood computed by importance sampling
+#>       -2LL= 614.4911 
+#>       AIC = 624.4911 
+#>       BIC = 622.5382 
+#> ----------------------------------------------------
f_saemix_fomc <- saemix(m_saemix_fomc, d_saemix_parent, saemix_options)
+
#> Running main SAEM algorithm
+#> [1] "Thu Nov  5 23:53:33 2020"
+#> ..
+#>     Minimisation finished
+#> [1] "Thu Nov  5 23:53:34 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> -----------------------------------
+#> ----          Data             ----
+#> -----------------------------------
+#> Object of class SaemixData
+#>     longitudinal data for use with the SAEM algorithm
+#> Dataset ds_saemix 
+#>     Structured data: value ~ time + name | ds 
+#>     X variable for graphs: time () 
+#> Dataset characteristics:
+#>     number of subjects:     5 
+#>     number of observations: 90 
+#>     average/min/max nb obs: 18.00  /  16  /  20 
+#> First 10 lines of data:
+#>           ds time   name value mdv cens occ ytype
+#> 1  Dataset 6    0 parent  97.2   0    0   1     1
+#> 2  Dataset 6    0 parent  96.4   0    0   1     1
+#> 3  Dataset 6    3 parent  71.1   0    0   1     1
+#> 4  Dataset 6    3 parent  69.2   0    0   1     1
+#> 5  Dataset 6    6 parent  58.1   0    0   1     1
+#> 6  Dataset 6    6 parent  56.6   0    0   1     1
+#> 7  Dataset 6   10 parent  44.4   0    0   1     1
+#> 8  Dataset 6   10 parent  43.4   0    0   1     1
+#> 9  Dataset 6   20 parent  33.3   0    0   1     1
+#> 10 Dataset 6   20 parent  29.2   0    0   1     1
+#> -----------------------------------
+#> ----          Model            ----
+#> -----------------------------------
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     psi[id, 1]/(xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 
+#>         2])
+#> }
+#> <bytecode: 0x55555998dc50>
+#> <environment: 0x5555595d7668>
+#>   Nb of parameters: 3 
+#>       parameter names:  parent_0 log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_alpha normal       Estimated
+#> [3,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           parent_0 log_alpha log_beta
+#> parent_0         1         0        0
+#> log_alpha        0         1        0
+#> log_beta         0         0        1
+#>   Error model: constant , initial values: a.1=1.91976382242696 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0  log_alpha log_beta
+#> Pop.CondInit 94.43884 -0.2226095 2.048192
+#> -----------------------------------
+#> ----    Key algorithm options  ----
+#> -----------------------------------
+#>     Estimation of individual parameters (MAP)
+#>     Estimation of standard errors and linearised log-likelihood
+#>     Estimation of log-likelihood by importance sampling
+#>     Number of iterations:  K1=200, K2=80 
+#>     Number of chains:  10 
+#>     Seed:  123456 
+#>     Number of MCMC iterations for IS:  5000 
+#>     Simulations:
+#>         nb of simulated datasets used for npde:  1000 
+#>         nb of simulated datasets used for VPC:  100 
+#>     Input/output
+#>         save the results to a file:  FALSE 
+#>         save the graphs to files:  FALSE 
+#> ----------------------------------------------------
+#> ----                  Results                   ----
+#> ----------------------------------------------------
+#> -----------------  Fixed effects  ------------------
+#> ----------------------------------------------------
+#>   Parameter Estimate SE   CV(%)
+#>   parent_0  94.49    1.18   1.2
+#>   log_alpha -0.21    0.30 142.0
+#>   log_beta   2.06    0.21  10.4
+#> a a.1        2.28    0.19   8.2
+#> ----------------------------------------------------
+#> -----------  Variance of random effects  -----------
+#> ----------------------------------------------------
+#>           Parameter        Estimate SE   CV(%)
+#> parent_0  omega2.parent_0  4.66     4.34 93   
+#> log_alpha omega2.log_alpha 0.45     0.29 65   
+#> log_beta  omega2.log_beta  0.19     0.14 75   
+#> ----------------------------------------------------
+#> ------  Correlation matrix of random effects  ------
+#> ----------------------------------------------------
+#>                  omega2.parent_0 omega2.log_alpha omega2.log_beta
+#> omega2.parent_0  1               0                0              
+#> omega2.log_alpha 0               1                0              
+#> omega2.log_beta  0               0                1              
+#> ----------------------------------------------------
+#> ---------------  Statistical criteria  -------------
+#> ----------------------------------------------------
+#> Likelihood computed by linearisation
+#>       -2LL= 454.0598 
+#>       AIC = 468.0598 
+#>       BIC = 465.3259 
+#> 
+#> Likelihood computed by importance sampling
+#>       -2LL= 453.7499 
+#>       AIC = 467.7499 
+#>       BIC = 465.016 
+#> ----------------------------------------------------
f_saemix_dfop <- saemix(m_saemix_dfop, d_saemix_parent, saemix_options)
+
#> Running main SAEM algorithm
+#> [1] "Thu Nov  5 23:53:35 2020"
+#> ..
+#>     Minimisation finished
+#> [1] "Thu Nov  5 23:53:37 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> -----------------------------------
+#> ----          Data             ----
+#> -----------------------------------
+#> Object of class SaemixData
+#>     longitudinal data for use with the SAEM algorithm
+#> Dataset ds_saemix 
+#>     Structured data: value ~ time + name | ds 
+#>     X variable for graphs: time () 
+#> Dataset characteristics:
+#>     number of subjects:     5 
+#>     number of observations: 90 
+#>     average/min/max nb obs: 18.00  /  16  /  20 
+#> First 10 lines of data:
+#>           ds time   name value mdv cens occ ytype
+#> 1  Dataset 6    0 parent  97.2   0    0   1     1
+#> 2  Dataset 6    0 parent  96.4   0    0   1     1
+#> 3  Dataset 6    3 parent  71.1   0    0   1     1
+#> 4  Dataset 6    3 parent  69.2   0    0   1     1
+#> 5  Dataset 6    6 parent  58.1   0    0   1     1
+#> 6  Dataset 6    6 parent  56.6   0    0   1     1
+#> 7  Dataset 6   10 parent  44.4   0    0   1     1
+#> 8  Dataset 6   10 parent  43.4   0    0   1     1
+#> 9  Dataset 6   20 parent  33.3   0    0   1     1
+#> 10 Dataset 6   20 parent  29.2   0    0   1     1
+#> -----------------------------------
+#> ----          Model            ----
+#> -----------------------------------
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     g <- plogis(psi[id, 4])
+#>     t = xidep[, "time"]
+#>     psi[id, 1] * (g * exp(-exp(psi[id, 2]) * t) + (1 - g) * exp(-exp(psi[id, 
+#>         3]) * t))
+#> }
+#> <bytecode: 0x55555998e548>
+#> <environment: 0x555558225bf0>
 #>   Nb of parameters: 4 
-#>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
+#>       parameter names:  parent_0 log_k1 log_k2 g_qlogis 
 #>       distribution:
-#>      Parameter      Distribution Estimated
-#> [1,] parent_0       normal       Estimated
-#> [2,] log_k_parent   normal       Estimated
-#> [3,] log_k_A1       normal       Estimated
-#> [4,] f_parent_ilr_1 normal       Estimated
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_k1    normal       Estimated
+#> [3,] log_k2    normal       Estimated
+#> [4,] g_qlogis  normal       Estimated
 #>   Variance-covariance matrix:
-#>                parent_0 log_k_parent log_k_A1 f_parent_ilr_1
-#> parent_0              1            0        0              0
-#> log_k_parent          0            1        0              0
-#> log_k_A1              0            0        1              0
-#> f_parent_ilr_1        0            0        0              1
-#>   Error model: constant , initial values: a.1=4.97259024646577 
+#>          parent_0 log_k1 log_k2 g_qlogis
+#> parent_0        1      0      0        0
+#> log_k1          0      1      0        0
+#> log_k2          0      0      1        0
+#> g_qlogis        0      0      0        1
+#>   Error model: constant , initial values: a.1=1.94671278396371 
 #>     No covariate in the model.
 #>     Initial values
-#>              parent_0 log_k_parent  log_k_A1 f_parent_ilr_1
-#> Pop.CondInit 86.53449    -3.207005 -3.060308      -1.920449
+#>              parent_0    log_k1    log_k2 g_qlogis
+#> Pop.CondInit 94.08322 -1.834163 -4.210797  0.11002
 #> -----------------------------------
 #> ----    Key algorithm options  ----
 #> -----------------------------------
@@ -383,69 +745,195 @@ variances of the deviations of the parameters from these mean values.

 #> ----------------------------------------------------
 #> -----------------  Fixed effects  ------------------
 #> ----------------------------------------------------
-#>   Parameter      Estimate SE   CV(%)
-#>   parent_0       86.09    1.57  1.8 
-#>   log_k_parent   -3.21    0.59 18.5 
-#>   log_k_A1       -4.69    0.31  6.6 
-#>   f_parent_ilr_1 -0.34    0.30 89.2 
-#> a a.1             4.69    0.27  5.8 
+#>   Parameter Estimate SE   CV(%)
+#>   parent_0  93.97    1.35   1.4
+#>   log_k1    -2.37    0.58  24.5
+#>   log_k2    -3.63    0.87  24.0
+#>   g_qlogis  -0.14    0.34 246.1
+#> a a.1        2.32    0.19   8.3
 #> ----------------------------------------------------
 #> -----------  Variance of random effects  -----------
 #> ----------------------------------------------------
-#>                Parameter             Estimate SE   CV(%)
-#> parent_0       omega2.parent_0       7.07     7.72 109  
-#> log_k_parent   omega2.log_k_parent   1.75     1.11  63  
-#> log_k_A1       omega2.log_k_A1       0.28     0.28  99  
-#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.39     0.27  71  
+#>          Parameter       Estimate SE   CV(%)
+#> parent_0 omega2.parent_0 6.97     5.72  82  
+#> log_k1   omega2.log_k1   1.63     1.06  65  
+#> log_k2   omega2.log_k2   3.73     2.39  64  
+#> g_qlogis omega2.g_qlogis 0.16     0.27 173  
 #> ----------------------------------------------------
 #> ------  Correlation matrix of random effects  ------
 #> ----------------------------------------------------
-#>                       omega2.parent_0 omega2.log_k_parent omega2.log_k_A1
-#> omega2.parent_0       1               0                   0              
-#> omega2.log_k_parent   0               1                   0              
-#> omega2.log_k_A1       0               0                   1              
-#> omega2.f_parent_ilr_1 0               0                   0              
-#>                       omega2.f_parent_ilr_1
-#> omega2.parent_0       0                    
-#> omega2.log_k_parent   0                    
-#> omega2.log_k_A1       0                    
-#> omega2.f_parent_ilr_1 1                    
+#>                 omega2.parent_0 omega2.log_k1 omega2.log_k2 omega2.g_qlogis
+#> omega2.parent_0 1               0             0             0              
+#> omega2.log_k1   0               1             0             0              
+#> omega2.log_k2   0               0             1             0              
+#> omega2.g_qlogis 0               0             0             1              
 #> ----------------------------------------------------
 #> ---------------  Statistical criteria  -------------
 #> ----------------------------------------------------
 #> Likelihood computed by linearisation
-#>       -2LL= 1064.35 
-#>       AIC = 1082.35 
-#>       BIC = 1078.835 
+#>       -2LL= 485.4627 
+#>       AIC = 503.4627 
+#>       BIC = 499.9477 
 #> 
 #> Likelihood computed by importance sampling
-#>       -2LL= 1063.475 
-#>       AIC = 1081.475 
-#>       BIC = 1077.96 
-#> ----------------------------------------------------
plot(f_saemix, plot.type = "convergence")
-
#> Plotting convergence plots
# }
-# Synthetic data with two-component error
-sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-dt50_sfo_in <- c(80, 90, 100, 111.111, 125)
-k_in <- log(2) / dt50_sfo_in
-
-SFO <- mkinmod(parent = mkinsub("SFO"))
-
-pred_sfo <- function(k) {
-  mkinpredict(SFO, c(k_parent = k),
-    c(parent = 100), sampling_times)
-}
-
-ds_sfo_mean <- lapply(k_in, pred_sfo)
-set.seed(123456L)
-ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) {
-  add_err(ds, sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
-     n = 1)[[1]]
- })
-# \dontrun{
-f_mmkin_syn <- mmkin("SFO", ds_sfo_syn, error_model = "tc", quiet = TRUE)
-# plot(f_mmkin_syn)
-m_saemix_tc <- saemix_model(f_mmkin_syn, cores = 1)
+#>       -2LL= 473.563 
+#>       AIC = 491.563 
+#>       BIC = 488.048 
+#> ----------------------------------------------------
compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
+
#> Likelihoods computed by importance sampling 
#>        AIC      BIC
+#> 1 624.4911 622.5382
+#> 2 467.7499 465.0160
+#> 3 491.5630 488.0480
f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
+m_saemix_fomc_tc <- saemix_model(f_mmkin_parent_tc["FOMC", ])
+
#> 
+#> 
+#> The following SaemixModel object was successfully created:
+#> 
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     psi[id, 1]/(xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 
+#>         2])
+#> }
+#> <bytecode: 0x55555998dc50>
+#> <environment: 0x555559a957f8>
+#>   Nb of parameters: 3 
+#>       parameter names:  parent_0 log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_alpha normal       Estimated
+#> [3,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           parent_0 log_alpha log_beta
+#> parent_0         1         0        0
+#> log_alpha        0         1        0
+#> log_beta         0         0        1
+#>   Error model: combined , initial values: a.1=1.10728182011691 b.1=0.024889924291374 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0  log_alpha log_beta
+#> Pop.CondInit 93.13042 -0.1215336 2.230815
f_saemix_fomc_tc <- saemix(m_saemix_fomc_tc, d_saemix_parent, saemix_options)
+
#> Running main SAEM algorithm
+#> [1] "Thu Nov  5 23:53:38 2020"
+#> ..
+#>     Minimisation finished
+#> [1] "Thu Nov  5 23:53:42 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> -----------------------------------
+#> ----          Data             ----
+#> -----------------------------------
+#> Object of class SaemixData
+#>     longitudinal data for use with the SAEM algorithm
+#> Dataset ds_saemix 
+#>     Structured data: value ~ time + name | ds 
+#>     X variable for graphs: time () 
+#> Dataset characteristics:
+#>     number of subjects:     5 
+#>     number of observations: 90 
+#>     average/min/max nb obs: 18.00  /  16  /  20 
+#> First 10 lines of data:
+#>           ds time   name value mdv cens occ ytype
+#> 1  Dataset 6    0 parent  97.2   0    0   1     1
+#> 2  Dataset 6    0 parent  96.4   0    0   1     1
+#> 3  Dataset 6    3 parent  71.1   0    0   1     1
+#> 4  Dataset 6    3 parent  69.2   0    0   1     1
+#> 5  Dataset 6    6 parent  58.1   0    0   1     1
+#> 6  Dataset 6    6 parent  56.6   0    0   1     1
+#> 7  Dataset 6   10 parent  44.4   0    0   1     1
+#> 8  Dataset 6   10 parent  43.4   0    0   1     1
+#> 9  Dataset 6   20 parent  33.3   0    0   1     1
+#> 10 Dataset 6   20 parent  29.2   0    0   1     1
+#> -----------------------------------
+#> ----          Model            ----
+#> -----------------------------------
+#> Nonlinear mixed-effects model
+#>   Model function:  Mixed model generated from mmkin object  Model type:  structural
+#> function (psi, id, xidep) 
+#> {
+#>     psi[id, 1]/(xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 
+#>         2])
+#> }
+#> <bytecode: 0x55555998dc50>
+#> <environment: 0x555559a957f8>
+#>   Nb of parameters: 3 
+#>       parameter names:  parent_0 log_alpha log_beta 
+#>       distribution:
+#>      Parameter Distribution Estimated
+#> [1,] parent_0  normal       Estimated
+#> [2,] log_alpha normal       Estimated
+#> [3,] log_beta  normal       Estimated
+#>   Variance-covariance matrix:
+#>           parent_0 log_alpha log_beta
+#> parent_0         1         0        0
+#> log_alpha        0         1        0
+#> log_beta         0         0        1
+#>   Error model: combined , initial values: a.1=1.10728182011691 b.1=0.024889924291374 
+#>     No covariate in the model.
+#>     Initial values
+#>              parent_0  log_alpha log_beta
+#> Pop.CondInit 93.13042 -0.1215336 2.230815
+#> -----------------------------------
+#> ----    Key algorithm options  ----
+#> -----------------------------------
+#>     Estimation of individual parameters (MAP)
+#>     Estimation of standard errors and linearised log-likelihood
+#>     Estimation of log-likelihood by importance sampling
+#>     Number of iterations:  K1=200, K2=80 
+#>     Number of chains:  10 
+#>     Seed:  123456 
+#>     Number of MCMC iterations for IS:  5000 
+#>     Simulations:
+#>         nb of simulated datasets used for npde:  1000 
+#>         nb of simulated datasets used for VPC:  100 
+#>     Input/output
+#>         save the results to a file:  FALSE 
+#>         save the graphs to files:  FALSE 
+#> ----------------------------------------------------
+#> ----                  Results                   ----
+#> ----------------------------------------------------
+#> -----------------  Fixed effects  ------------------
+#> ----------------------------------------------------
+#>   Parameter Estimate SE     CV(%)
+#>   parent_0  94.4481  1.2052   1.3
+#>   log_alpha -0.2088  0.3059 146.5
+#>   log_beta   2.0668  0.2182  10.6
+#> a a.1        2.4273  0.3178  13.1
+#> b b.1       -0.0037  0.0062 168.3
+#> ----------------------------------------------------
+#> -----------  Variance of random effects  -----------
+#> ----------------------------------------------------
+#>           Parameter        Estimate SE   CV(%)
+#> parent_0  omega2.parent_0  5.34     4.58 86   
+#> log_alpha omega2.log_alpha 0.46     0.29 65   
+#> log_beta  omega2.log_beta  0.20     0.15 74   
+#> ----------------------------------------------------
+#> ------  Correlation matrix of random effects  ------
+#> ----------------------------------------------------
+#>                  omega2.parent_0 omega2.log_alpha omega2.log_beta
+#> omega2.parent_0  1               0                0              
+#> omega2.log_alpha 0               1                0              
+#> omega2.log_beta  0               0                1              
+#> ----------------------------------------------------
+#> ---------------  Statistical criteria  -------------
+#> ----------------------------------------------------
+#> Likelihood computed by linearisation
+#>       -2LL= 453.7703 
+#>       AIC = 469.7703 
+#>       BIC = 466.6458 
+#> 
+#> Likelihood computed by importance sampling
+#>       -2LL= 453.6186 
+#>       AIC = 469.6186 
+#>       BIC = 466.4942 
+#> ----------------------------------------------------
compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
+
#> Likelihoods computed by importance sampling 
#>        AIC      BIC
+#> 1 467.7499 465.0160
+#> 2 469.6186 466.4942

+dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+  A1 = mkinsub("SFO"))
+
#> Successfully compiled differential equation model from auto-generated C code.
f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
+m_saemix <- saemix_model(f_mmkin)
 
#> 
 #> 
 #> The following SaemixModel object was successfully created:
@@ -468,7 +956,7 @@ variances of the deviations of the parameters from these mean values.

 #>             transform_fractions = object[[1]]$transform_fractions)
 #>         odeparms <- c(odeparms_optim, odeparms_fixed)
 #>         xidep_i <- subset(xidep, id == i)
-#>         if (analytical) {
+#>         if (solution_type == "analytical") {
 #>             out_values <- mkin_model$deg_func(xidep_i, odeini, 
 #>                 odeparms)
 #>         }
@@ -476,7 +964,7 @@ variances of the deviations of the parameters from these mean values.

 #>             i_time <- xidep_i$time
 #>             i_name <- xidep_i$name
 #>             out_wide <- mkinpredict(mkin_model, odeparms = odeparms, 
-#>                 odeini = odeini, solution_type = object[[1]]$solution_type, 
+#>                 odeini = odeini, solution_type = solution_type, 
 #>                 outtimes = sort(unique(i_time)))
 #>             out_index <- cbind(as.character(i_time), as.character(i_name))
 #>             out_values <- out_wide[out_index]
@@ -486,23 +974,31 @@ variances of the deviations of the parameters from these mean values.

 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x55555d62aeb8>
-#> <environment: 0x55555cd8e028>
-#>   Nb of parameters: 2 
-#>       parameter names:  parent_0 log_k_parent 
+#> <bytecode: 0x55555998cba0>
+#> <environment: 0x55555bd1fee8>
+#>   Nb of parameters: 6 
+#>       parameter names:  parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 g_qlogis 
 #>       distribution:
-#>      Parameter    Distribution Estimated
-#> [1,] parent_0     normal       Estimated
-#> [2,] log_k_parent normal       Estimated
+#>      Parameter       Distribution Estimated
+#> [1,] parent_0        normal       Estimated
+#> [2,] log_k_A1        normal       Estimated
+#> [3,] f_parent_qlogis normal       Estimated
+#> [4,] log_k1          normal       Estimated
+#> [5,] log_k2          normal       Estimated
+#> [6,] g_qlogis        normal       Estimated
 #>   Variance-covariance matrix:
-#>              parent_0 log_k_parent
-#> parent_0            1            0
-#> log_k_parent        0            1
-#>   Error model: combined , initial values: a.1=1.05209877924905 b.1=0.0586479225303944 
+#>                 parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 g_qlogis
+#> parent_0               1        0               0      0      0        0
+#> log_k_A1               0        1               0      0      0        0
+#> f_parent_qlogis        0        0               1      0      0        0
+#> log_k1                 0        0               0      1      0        0
+#> log_k2                 0        0               0      0      1        0
+#> g_qlogis               0        0               0      0      0        1
+#>   Error model: constant , initial values: a.1=1.64723790168612 
 #>     No covariate in the model.
 #>     Initial values
-#>              parent_0 log_k_parent
-#> Pop.CondInit  100.315    -4.962075
d_saemix_tc <- saemix_data(f_mmkin_syn)
+#>              parent_0  log_k_A1 f_parent_qlogis    log_k1    log_k2  g_qlogis
+#> Pop.CondInit 93.81015 -9.764746      -0.9711148 -1.879937 -4.270814 0.1356441
d_saemix <- saemix_data(f_mmkin)
 
#> 
 #> 
 #> The following SaemixData object was successfully created:
@@ -511,12 +1007,12 @@ variances of the deviations of the parameters from these mean values.

 #>     longitudinal data for use with the SAEM algorithm
 #> Dataset ds_saemix 
 #>     Structured data: value ~ time + name | ds 
-#>     X variable for graphs: time () 
f_saemix_tc <- saemix(m_saemix_tc, d_saemix_tc, saemix_options)
+#>     X variable for graphs: time () 
f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
 
#> Running main SAEM algorithm
-#> [1] "Thu Nov  5 08:28:50 2020"
+#> [1] "Thu Nov  5 23:53:43 2020"
 #> ..
 #>     Minimisation finished
-#> [1] "Thu Nov  5 08:29:41 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> [1] "Thu Nov  5 23:56:33 2020"
#> Nonlinear mixed-effects model fit by the SAEM algorithm
 #> -----------------------------------
 #> ----          Data             ----
 #> -----------------------------------
@@ -527,20 +1023,20 @@ variances of the deviations of the parameters from these mean values.

 #>     X variable for graphs: time () 
 #> Dataset characteristics:
 #>     number of subjects:     5 
-#>     number of observations: 90 
-#>     average/min/max nb obs: 18.00  /  18  /  18 
+#>     number of observations: 170 
+#>     average/min/max nb obs: 34.00  /  30  /  38 
 #> First 10 lines of data:
-#>    ds time   name value mdv cens occ ytype
-#> 1   1    0 parent 105.9   0    0   1     1
-#> 2   1    0 parent  98.0   0    0   1     1
-#> 3   1    1 parent  96.6   0    0   1     1
-#> 4   1    1 parent  99.8   0    0   1     1
-#> 5   1    3 parent 113.0   0    0   1     1
-#> 6   1    3 parent 103.2   0    0   1     1
-#> 7   1    7 parent 102.9   0    0   1     1
-#> 8   1    7 parent 110.8   0    0   1     1
-#> 9   1   14 parent  95.9   0    0   1     1
-#> 10  1   14 parent  85.9   0    0   1     1
+#>           ds time   name value mdv cens occ ytype
+#> 1  Dataset 6    0 parent  97.2   0    0   1     1
+#> 2  Dataset 6    0 parent  96.4   0    0   1     1
+#> 3  Dataset 6    3 parent  71.1   0    0   1     1
+#> 4  Dataset 6    3 parent  69.2   0    0   1     1
+#> 5  Dataset 6    6 parent  58.1   0    0   1     1
+#> 6  Dataset 6    6 parent  56.6   0    0   1     1
+#> 7  Dataset 6   10 parent  44.4   0    0   1     1
+#> 8  Dataset 6   10 parent  43.4   0    0   1     1
+#> 9  Dataset 6   20 parent  33.3   0    0   1     1
+#> 10 Dataset 6   20 parent  29.2   0    0   1     1
 #> -----------------------------------
 #> ----          Model            ----
 #> -----------------------------------
@@ -562,7 +1058,7 @@ variances of the deviations of the parameters from these mean values.

 #>             transform_fractions = object[[1]]$transform_fractions)
 #>         odeparms <- c(odeparms_optim, odeparms_fixed)
 #>         xidep_i <- subset(xidep, id == i)
-#>         if (analytical) {
+#>         if (solution_type == "analytical") {
 #>             out_values <- mkin_model$deg_func(xidep_i, odeini, 
 #>                 odeparms)
 #>         }
@@ -570,7 +1066,7 @@ variances of the deviations of the parameters from these mean values.

 #>             i_time <- xidep_i$time
 #>             i_name <- xidep_i$name
 #>             out_wide <- mkinpredict(mkin_model, odeparms = odeparms, 
-#>                 odeini = odeini, solution_type = object[[1]]$solution_type, 
+#>                 odeini = odeini, solution_type = solution_type, 
 #>                 outtimes = sort(unique(i_time)))
 #>             out_index <- cbind(as.character(i_time), as.character(i_name))
 #>             out_values <- out_wide[out_index]
@@ -580,23 +1076,31 @@ variances of the deviations of the parameters from these mean values.

 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x55555d62aeb8>
-#> <environment: 0x55555cd8e028>
-#>   Nb of parameters: 2 
-#>       parameter names:  parent_0 log_k_parent 
+#> <bytecode: 0x55555998cba0>
+#> <environment: 0x55555bd1fee8>
+#>   Nb of parameters: 6 
+#>       parameter names:  parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 g_qlogis 
 #>       distribution:
-#>      Parameter    Distribution Estimated
-#> [1,] parent_0     normal       Estimated
-#> [2,] log_k_parent normal       Estimated
+#>      Parameter       Distribution Estimated
+#> [1,] parent_0        normal       Estimated
+#> [2,] log_k_A1        normal       Estimated
+#> [3,] f_parent_qlogis normal       Estimated
+#> [4,] log_k1          normal       Estimated
+#> [5,] log_k2          normal       Estimated
+#> [6,] g_qlogis        normal       Estimated
 #>   Variance-covariance matrix:
-#>              parent_0 log_k_parent
-#> parent_0            1            0
-#> log_k_parent        0            1
-#>   Error model: combined , initial values: a.1=1.05209877924905 b.1=0.0586479225303944 
+#>                 parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 g_qlogis
+#> parent_0               1        0               0      0      0        0
+#> log_k_A1               0        1               0      0      0        0
+#> f_parent_qlogis        0        0               1      0      0        0
+#> log_k1                 0        0               0      1      0        0
+#> log_k2                 0        0               0      0      1        0
+#> g_qlogis               0        0               0      0      0        1
+#>   Error model: constant , initial values: a.1=1.64723790168612 
 #>     No covariate in the model.
 #>     Initial values
-#>              parent_0 log_k_parent
-#> Pop.CondInit  100.315    -4.962075
+#>              parent_0  log_k_A1 f_parent_qlogis    log_k1    log_k2  g_qlogis
+#> Pop.CondInit 93.81015 -9.764746      -0.9711148 -1.879937 -4.270814 0.1356441
 #> -----------------------------------
 #> ----    Key algorithm options  ----
 #> -----------------------------------
@@ -618,37 +1122,55 @@ variances of the deviations of the parameters from these mean values.

 #> ----------------------------------------------------
 #> -----------------  Fixed effects  ------------------
 #> ----------------------------------------------------
-#>   Parameter    Estimate SE    CV(%) 
-#>   parent_0     100.232  1.266    1.3
-#>   log_k_parent  -4.961  0.089    1.8
-#> a a.1           -0.106  1.211 1142.0
-#> b b.1            0.071  0.017   24.2
+#>   Parameter       Estimate SE   CV(%)
+#>   parent_0        93.78    1.35   1.4
+#>   log_k_A1        -6.05    1.12  18.5
+#>   f_parent_qlogis -0.97    0.20  21.1
+#>   log_k1          -2.46    0.51  20.7
+#>   log_k2          -3.63    0.95  26.3
+#>   g_qlogis        -0.08    0.36 447.7
+#> a a.1              1.88    0.11   5.9
 #> ----------------------------------------------------
 #> -----------  Variance of random effects  -----------
 #> ----------------------------------------------------
-#>              Parameter           Estimate SE    CV(%)
-#> parent_0     omega2.parent_0     3.334    5.024 151  
-#> log_k_parent omega2.log_k_parent 0.036    0.024  68  
+#>                 Parameter              Estimate SE   CV(%)
+#> parent_0        omega2.parent_0        7.85     5.76  73  
+#> log_k_A1        omega2.log_k_A1        4.27     3.44  80  
+#> f_parent_qlogis omega2.f_parent_qlogis 0.20     0.13  65  
+#> log_k1          omega2.log_k1          1.08     0.77  72  
+#> log_k2          omega2.log_k2          4.24     2.83  67  
+#> g_qlogis        omega2.g_qlogis        0.21     0.26 123  
 #> ----------------------------------------------------
 #> ------  Correlation matrix of random effects  ------
 #> ----------------------------------------------------
-#>                     omega2.parent_0 omega2.log_k_parent
-#> omega2.parent_0     1               0                  
-#> omega2.log_k_parent 0               1                  
+#>                        omega2.parent_0 omega2.log_k_A1 omega2.f_parent_qlogis
+#> omega2.parent_0        1               0               0                     
+#> omega2.log_k_A1        0               1               0                     
+#> omega2.f_parent_qlogis 0               0               1                     
+#> omega2.log_k1          0               0               0                     
+#> omega2.log_k2          0               0               0                     
+#> omega2.g_qlogis        0               0               0                     
+#>                        omega2.log_k1 omega2.log_k2 omega2.g_qlogis
+#> omega2.parent_0        0             0             0              
+#> omega2.log_k_A1        0             0             0              
+#> omega2.f_parent_qlogis 0             0             0              
+#> omega2.log_k1          1             0             0              
+#> omega2.log_k2          0             1             0              
+#> omega2.g_qlogis        0             0             1              
 #> ----------------------------------------------------
 #> ---------------  Statistical criteria  -------------
 #> ----------------------------------------------------
 #> Likelihood computed by linearisation
-#>       -2LL= 575.5586 
-#>       AIC = 587.5586 
-#>       BIC = 585.2153 
+#>       -2LL= 879.7721 
+#>       AIC = 905.7721 
+#>       BIC = 900.6948 
 #> 
 #> Likelihood computed by importance sampling
-#>       -2LL= 575.7797 
-#>       AIC = 587.7797 
-#>       BIC = 585.4364 
-#> ----------------------------------------------------
plot(f_saemix_tc, plot.type = "convergence")
-
#> Plotting convergence plots
# }
+#>       -2LL= 816.8276 
+#>       AIC = 842.8276 
+#>       BIC = 837.7503 
+#> ----------------------------------------------------

+# }