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: 0x55555e4213d8>
+#> <environment: 0x55555c47cdf0>
+#>   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] "Wed May 27 05:55:50 2020"
+#> ..
+#>     Minimisation finished
+#> [1] "Wed May 27 06:01:54 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: 0x55555e4213d8>
+#> <environment: 0x55555c47cdf0>
+#>   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
# }
+

+ +