Create saemix models from mmkin row objects

- - -

This function sets up a nonlinear mixed effects model for an mmkin row -object for use with the saemix package. An mmkin row object is essentially a -list of mkinfit objects that have been obtained by fitting the same model to -a list of datasets.

- -

saemix_model(object, cores = parallel::detectCores())
-
-saemix_data(object, ...)

- -

Arguments

- - - - - - - - - - - - - - -

object	An mmkin row object containing several fits of the same model to different datasets
cores	The number of cores to be used for multicore processing. -On Windows machines, cores > 1 is currently not supported.
...	Further parameters passed to saemix::saemixData

- -

Value

- -

An saemix::SaemixModel object.

An saemix::SaemixData object.

Details

- -

Starting values for the fixed effects (population mean parameters, argument psi0 of -saemix::saemixModel() are the mean values of the parameters found using -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)
-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)
#> 
-#> 
-#> 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) 
-#> {
-#>     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)
-#> }
-#> <bytecode: 0x555559668108>
-#> <environment: 0x555559677c08>
-#>   Nb of parameters: 4 
-#>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
-#>       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
-#>   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=1 
-#>     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)
#> 
-#> 
-#> 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,
-  save = FALSE, save.graphs = FALSE, displayProgress = FALSE,
-  nbiter.saemix = c(200, 80))
-f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
#> Running main SAEM algorithm
-#> [1] "Tue May 26 18:25:16 2020"
-#> ..
-#>     Minimisation finished
-#> [1] "Tue May 26 18:31:09 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: 170 
-#>     average/min/max nb obs: 34.00  /  30  /  38 
-#> 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) 
-#> {
-#>     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)
-#> }
-#> <bytecode: 0x555559668108>
-#> <environment: 0x555559677c08>
-#>   Nb of parameters: 4 
-#>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
-#>       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
-#>   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=1 
-#>     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
-#> -----------------------------------
-#> ----    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(%)
-#> [1,] parent_0       86.14    1.61  1.9 
-#> [2,] log_k_parent   -3.21    0.59 18.5 
-#> [3,] log_k_A1       -4.66    0.30  6.4 
-#> [4,] f_parent_ilr_1 -0.33    0.30 91.7 
-#> [5,] a.1             4.68    0.27  5.8 
-#> ----------------------------------------------------
-#> -----------  Variance of random effects  -----------
-#> ----------------------------------------------------
-#>                Parameter             Estimate SE   CV(%)
-#> parent_0       omega2.parent_0       7.71     8.14 106  
-#> log_k_parent   omega2.log_k_parent   1.76     1.12  63  
-#> log_k_A1       omega2.log_k_A1       0.26     0.26 101  
-#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.39     0.28  71  
-#> ----------------------------------------------------
-#> ------  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                    
-#> ----------------------------------------------------
-#> ---------------  Statistical criteria  -------------
-#> ----------------------------------------------------
-#> Likelihood computed by linearisation
-#>       -2LL= 1064.364 
-#>       AIC = 1082.364 
-#>       BIC = 1078.848 
-#> 
-#> Likelihood computed by importance sampling
-#>       -2LL= 1063.462 
-#>       AIC = 1081.462 
-#>       BIC = 1077.947 
-#> ----------------------------------------------------
plot(f_saemix, plot.type = "convergence")
#> Plotting convergence plots
# }
-

- -