From 675a733fa2acc08daabb9b8b571c7d658f281f73 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 26 May 2020 18:38:51 +0200 Subject: Use all cores per default, confint tolerance Also, use more intelligent starting values for the variance of the random effects for saemix. While this does not appear to speed up the convergence, it shows where this variance is greatly reduced by using mixed-effects models as opposed to the separate independent fits. --- docs/reference/confint.mkinfit.html | 20 ++++-- docs/reference/mmkin.html | 16 +++-- docs/reference/nlme.html | 4 +- docs/reference/parms.html | 6 +- docs/reference/saemix-1.png | Bin 37443 -> 31551 bytes docs/reference/saemix-2.png | Bin 38557 -> 58815 bytes docs/reference/saemix.html | 135 ++++++++++++------------------------ 7 files changed, 77 insertions(+), 104 deletions(-) (limited to 'docs/reference') diff --git a/docs/reference/confint.mkinfit.html b/docs/reference/confint.mkinfit.html index 190494bc..0686c7bb 100644 --- a/docs/reference/confint.mkinfit.html +++ b/docs/reference/confint.mkinfit.html @@ -79,7 +79,7 @@ method of Venzon and Moolgavkar (1988)." /> mkin - 0.9.50.2 + 0.9.50.3 @@ -116,6 +116,9 @@ method of Venzon and Moolgavkar (1988)." />
  • Example evaluation of NAFTA SOP Attachment examples
    • +
  • + Some benchmark timings +
  • @@ -168,7 +171,8 @@ method of Venzon and Moolgavkar (1988).

    method = c("quadratic", "profile"), transformed = TRUE, backtransform = TRUE, - cores = round(detectCores()/2), + cores = parallel::detectCores(), + rel_tol = 0.01, quiet = FALSE, ... ) @@ -222,6 +226,12 @@ their confidence intervals?

    cores

    The number of cores to be used for multicore processing. On Windows machines, cores > 1 is currently not supported.

    + + + rel_tol +

    If the method is 'profile', what should be the accuracy +of the lower and upper bounds, relative to the estimate obtained from +the quadratic method?

    quiet @@ -270,13 +280,13 @@ Profile-Likelihood Based Confidence Intervals, Applied Statistics, 37, SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
    #> User System verstrichen -#> 3.410 0.000 3.412
    # Using more cores does not save much time here, as parent_0 takes up most of the time +system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
    #> user system elapsed +#> 3.689 0.991 3.361
    # Using more cores does not save much time here, as parent_0 takes up most of the time # If we additionally exclude parent_0 (the confidence of which is often of # minor interest), we get a nice performance improvement from about 50 # seconds to about 12 seconds if we use at least four cores system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile", - c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
    #> Profiling the likelihood
    #> Warning: scheduled cores 1, 2, 3 encountered errors in user code, all values of the jobs will be affected
    #> Error in dimnames(x) <- dn: Länge von 'dimnames' [2] ungleich der Arrayausdehnung
    #> Timing stopped at: 0.008 0.044 0.201
    ci_profile
    #> 2.5% 97.5% + c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
    #> Profiling the likelihood
    #> Warning: scheduled cores 2, 1, 3 encountered errors in user code, all values of the jobs will be affected
    #> Error in dimnames(x) <- dn: length of 'dimnames' [2] not equal to array extent
    #> Timing stopped at: 0.007 0.042 0.193
    ci_profile
    #> 2.5% 97.5% #> parent_0 96.456003640 1.027703e+02 #> k_parent 0.090911032 1.071578e-01 #> k_m1 0.003892605 6.702778e-03 diff --git a/docs/reference/mmkin.html b/docs/reference/mmkin.html index 3be3b4b9..9628c017 100644 --- a/docs/reference/mmkin.html +++ b/docs/reference/mmkin.html @@ -75,7 +75,7 @@ datasets specified in its first two arguments." /> mkin - 0.9.50.2 + 0.9.50.3
    @@ -112,6 +112,9 @@ datasets specified in its first two arguments." />
  • Example evaluation of NAFTA SOP Attachment examples
    • +
  • + Some benchmark timings +
  • @@ -152,7 +155,7 @@ datasets specified in its first two arguments.

    mmkin(
   models = c("SFO", "FOMC", "DFOP"),
   datasets,
-  cores = round(detectCores()/2),
+  cores = detectCores(),
   cluster = NULL,
   ...
 )
    @@ -176,7 +179,8 @@ data for mkinfit.

    The number of cores to be used for multicore processing. This is only used when the cluster argument is NULL. On Windows machines, cores > 1 is not supported, you need to use the cluster -argument to use multiple logical processors.

    +argument to use multiple logical processors. Per default, all cores +detected by parallel::detectCores() are used.

    cluster @@ -215,9 +219,9 @@ plotting.

    time_default <- system.time(fits.0 <- mmkin(models, datasets, quiet = TRUE)) time_1 <- system.time(fits.4 <- mmkin(models, datasets, cores = 1, quiet = TRUE))
    #> Warning: Optimisation did not converge: #> false convergence (8)
    -time_default
    #> User System verstrichen -#> 4.520 0.374 1.284
    time_1
    #> User System verstrichen -#> 5.076 0.004 5.083
    +time_default
    #> user system elapsed +#> 4.457 0.561 1.328
    time_1
    #> user system elapsed +#> 5.031 0.004 5.038
    endpoints(fits.0[["SFO_lin", 2]])
    #> $ff #> parent_M1 parent_sink M1_M2 M1_sink #> 0.7340478 0.2659522 0.7505691 0.2494309 diff --git a/docs/reference/nlme.html b/docs/reference/nlme.html index b2d415dc..3462e52e 100644 --- a/docs/reference/nlme.html +++ b/docs/reference/nlme.html @@ -75,7 +75,7 @@ datasets. They are used internally by the nlme.mmkin() method." /> mkin - 0.9.50.2 + 0.9.50.3
    @@ -178,7 +178,7 @@ datasets. They are used internally by the nlme.m

    If random is FALSE (default), a named vector containing mean values of the fitted degradation model parameters. If random is TRUE, a list with fixed and random effects, in the format required by the start argument of -nlme for the case of a single grouping variable ds?

    +nlme for the case of a single grouping variable ds.

    A groupedData object

    See also

    diff --git a/docs/reference/parms.html b/docs/reference/parms.html index 432bbc88..2fe91c26 100644 --- a/docs/reference/parms.html +++ b/docs/reference/parms.html @@ -195,7 +195,8 @@ such matrices is returned.

    ds <- lapply(experimental_data_for_UBA_2019[6:10], function(x) subset(x$data[c("name", "time", "value")])) names(ds) <- paste("Dataset", 6:10) -fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) +# \dontrun{ +fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) parms(fits["SFO", ])
    #> Dataset 6 Dataset 7 Dataset 8 Dataset 9 Dataset 10 #> parent_0 88.52275400 82.666781678 86.8547308 91.7779306 82.14809450 #> k_parent_sink 0.05794659 0.009647805 0.2102974 0.1232258 0.00720421 @@ -259,7 +260,8 @@ such matrices is returned.

    #> log_k2 -3.5206791 -5.85402317 -2.5794240 -3.4233253 -5.676532 #> g_ilr -0.1463234 0.07627854 0.4719196 0.4477805 -0.460676 #> sigma 1.3569047 2.22130220 1.3416908 2.8715985 1.942068 -#>
    +#>
    # } +
    - 
    saemix_model(object)
+    
    saemix_model(object, cores = parallel::detectCores())
 
 saemix_data(object, ...)
    
 
@@ -164,6 +164,11 @@ a list of datasets.
    
       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
    
@@ -174,21 +179,22 @@ a list of datasets.
    
 
     
    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.
    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
    #> 
+  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:
 #> 
@@ -224,12 +230,12 @@ a list of datasets.
    
 #>             out_values <- out_wide[out_index]
 #>         }
 #>         return(out_values)
-#>     }, mc.cores = 15)
+#>     }, mc.cores = cores)
 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x555559875398>
-#> <environment: 0x55555973a248>
+#> <bytecode: 0x555559668108>
+#> <environment: 0x555559677c08>
 #>   Nb of parameters: 4 
 #>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
 #>       distribution:
@@ -248,8 +254,7 @@ a list of datasets.
    
 #>     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
-#> 
+#> Pop.CondInit 86.53449    -3.207005 -3.060308      -1.920449
    d_saemix <- saemix_data(f_mmkin)
    #> 
 #> 
 #> The following SaemixData object was successfully created:
 #> 
@@ -257,12 +262,14 @@ a list of datasets.
    
 #>     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"
-#> .
    #> .
    #> .
    #> .
    #> 
+#>     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] "Mon May 25 12:56:39 2020"
    #> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> [1] "Tue May 26 18:31:09 2020"
    #> Nonlinear mixed-effects model fit by the SAEM algorithm
 #> -----------------------------------
 #> ----          Data             ----
 #> -----------------------------------
@@ -322,12 +329,12 @@ a list of datasets.
    
 #>             out_values <- out_wide[out_index]
 #>         }
 #>         return(out_values)
-#>     }, mc.cores = 15)
+#>     }, mc.cores = cores)
 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x555559875398>
-#> <environment: 0x55555973a248>
+#> <bytecode: 0x555559668108>
+#> <environment: 0x555559677c08>
 #>   Nb of parameters: 4 
 #>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
 #>       distribution:
@@ -353,7 +360,7 @@ a list of datasets.
    
 #>     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 iterations:  K1=200, K2=80 
 #>     Number of chains:  10 
 #>     Seed:  123456 
 #>     Number of MCMC iterations for IS:  5000 
@@ -369,19 +376,19 @@ a list of datasets.
    
 #> -----------------  Fixed effects  ------------------
 #> ----------------------------------------------------
 #>      Parameter      Estimate SE   CV(%)
-#> [1,] parent_0       86.21    1.51  1.7 
+#> [1,] parent_0       86.14    1.61  1.9 
 #> [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 
+#> [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       6.07     7.08 117  
-#> log_k_parent   omega2.log_k_parent   1.75     1.11  63  
+#> 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.38     0.27  71  
+#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.39     0.28  71  
 #> ----------------------------------------------------
 #> ------  Correlation matrix of random effects  ------
 #> ----------------------------------------------------
@@ -399,66 +406,16 @@ a list of datasets.
    
 #> ---------------  Statistical criteria  -------------
 #> ----------------------------------------------------
 #> Likelihood computed by linearisation
-#>       -2LL= 1064.397 
-#>       AIC = 1082.397 
-#>       BIC = 1078.882 
+#>       -2LL= 1064.364 
+#>       AIC = 1082.364 
+#>       BIC = 1078.848 
 #> 
 #> 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 
    # }
-
    
+#>       -2LL= 1063.462 
+#>       AIC = 1081.462 
+#>       BIC = 1077.947 
+#> ----------------------------------------------------
    plot(f_saemix, plot.type = "convergence")
    #> Plotting convergence plots
    # }
+