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)
+
+saemix_data(object, ...)

+ +

Arguments

+ + + + + + + + + + +

object	An mmkin row object containing several fits of the same model to different datasets
...	Further parameters passed to saemix::saemixData

+ +

Value

+ +

An saemix::SaemixModel object.

An saemix::SaemixData object.

+ +

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.
f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE, cores = 5)
+# \dontrun{
+if (require(saemix)) {
+  m_saemix <- saemix_model(f_mmkin)
+  d_saemix <- saemix_data(f_mmkin)
+  saemix_options <- list(seed = 123456, save = FALSE, save.graphs = FALSE)
+    saemix(m_saemix, d_saemix, saemix_options)
+}
#> Loading required package: saemix
#> Package saemix, version 3.1.9000
+#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
#> 
+#> 
+#> 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 = 15)
+#>     res <- unlist(res_list)
+#>     return(res)
+#> }
+#> <bytecode: 0x555559875398>
+#> <environment: 0x55555973a248>
+#>   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
+#> 
+#> 
+#> 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 () 
+#> Running main SAEM algorithm
+#> [1] "Mon May 25 12:48:51 2020"
+#> .
#> .
#> .
#> .
#> 
+#>     Minimisation finished
+#> [1] "Mon May 25 12:56:39 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 = 15)
+#>     res <- unlist(res_list)
+#>     return(res)
+#> }
+#> <bytecode: 0x555559875398>
+#> <environment: 0x55555973a248>
+#>   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=300, K2=100 
+#>     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.21    1.51  1.7 
+#> [2,] log_k_parent   -3.21    0.59 18.5 
+#> [3,] log_k_A1       -4.64    0.29  6.3 
+#> [4,] f_parent_ilr_1 -0.32    0.30 93.2 
+#> [5,] a.1             4.69    0.27  5.8 
+#> ----------------------------------------------------
+#> -----------  Variance of random effects  -----------
+#> ----------------------------------------------------
+#>                Parameter             Estimate SE   CV(%)
+#> parent_0       omega2.parent_0       6.07     7.08 117  
+#> log_k_parent   omega2.log_k_parent   1.75     1.11  63  
+#> log_k_A1       omega2.log_k_A1       0.26     0.26 101  
+#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.38     0.27  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.397 
+#>       AIC = 1082.397 
+#>       BIC = 1078.882 
+#> 
+#> Likelihood computed by importance sampling
+#>       -2LL= 1063.161 
+#>       AIC = 1081.161 
+#>       BIC = 1077.646 
+#> ----------------------------------------------------
#> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> -----------------------------------------
+#> ----         Data and Model          ----
+#> -----------------------------------------
+#> Data
+#>     Dataset ds_saemix 
+#>     Longitudinal data: value ~ time + name | ds 
+#> 
+#> Model:
+#>     Mixed model generated from mmkin object 
+#>      4 parameters: parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
+#>      error model: constant 
+#>      No covariate
+#> 
+#> Key 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=300, K2=100 
+#>     Number of chains:  10 
+#>     Seed:  123456 
+#>     Number of MCMC iterations for IS:  5000 
+#>     Input/output
+#>         results not saved
+#>         no graphs
+#> ----------------------------------------------------
+#> ----                  Results                   ----
+#> Fixed effects
+#>  Parameter      Estimate   SE   CV(%)
+#>  parent_0       86.214   1.506  1.75 
+#>  log_k_parent   -3.210   0.593 18.47 
+#>  log_k_A1       -4.643   0.294  6.34 
+#>  f_parent_ilr_1 -0.322   0.300 93.24 
+#>  a.1             4.689   0.270  5.76 
+#> 
+#> Variance of random effects
+#>  Parameter             Estimate   SE   CV(%)
+#>  omega2.parent_0       6.068    7.078 116.7 
+#>  omega2.log_k_parent   1.752    1.111  63.4 
+#>  omega2.log_k_A1       0.256    0.257 100.5 
+#>  omega2.f_parent_ilr_1 0.385    0.273  70.8 
+#> 
+#> Statistical criteria
+#> Likelihood computed by linearisation
+#>       -2LL= 1064.397 
+#>        AIC= 1082.397 
+#>        BIC= 1078.882 
+#> Likelihood computed by importance sampling
+#>       -2LL= 1063.161 
+#>        AIC= 1081.161 
+#>        BIC= 1077.646 
# }
+

+ +