This function uses saemix::saemix()
as a backend for fitting nonlinear mixed
effects models created from mmkin row objects using the Stochastic Approximation
Expectation Maximisation algorithm (SAEM).
saem(object, control, ...) # S3 method for mmkin saem( object, control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE), cores = 1, verbose = FALSE, suppressPlot = TRUE, ... ) # S3 method for saem.mmkin print(x, digits = max(3, getOption("digits") - 3), ...) saemix_model(object, cores = 1, verbose = FALSE, ...) saemix_data(object, verbose = FALSE, ...)
object | An mmkin row object containing several fits of the same mkinmod model to different datasets |
---|---|
control | Passed to saemix::saemix |
... | Further parameters passed to saemix::saemixData and saemix::saemixModel. |
cores | The number of cores to be used for multicore processing using
|
verbose | Should we print information about created objects? |
suppressPlot | Should we suppress any plotting that is done by the saemix function? |
x | An saem.mmkin object to print |
digits | Number of digits to use for printing |
An S3 object of class 'saem.mmkin', containing the fitted saemix::SaemixObject as a list component named 'so'. The object also inherits from 'mixed.mmkin'.
An saemix::SaemixModel object.
An saemix::SaemixData object.
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 using mkinfit.
Starting values for the fixed effects (population mean parameters, argument
psi0 of saemix::saemixModel()
are the mean values of the parameters found
using mmkin.
# \dontrun{ 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) f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:21 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:23 2020"f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) f_saem_sfo <- saem(f_mmkin_parent["SFO", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:24 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:26 2020"f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:26 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:28 2020"f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:29 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:32 2020"# The returned saem.mmkin object contains an SaemixObject, therefore we can use # functions from saemix library(saemix)#>#>#> Likelihoods computed by importance sampling#> AIC BIC #> 1 624.2428 622.2900 #> 2 467.7644 465.0305 #> 3 491.3541 487.8391#> Plotting convergence plots#> Plotting individual fits#> Simulating data using nsim = 1000 simulated datasets #> Computing WRES and npde . #> Plotting npde#> --------------------------------------------- #> Distribution of npde: #> mean= -0.01736 (SE= 0.098 ) #> variance= 0.8562 (SE= 0.13 ) #> skewness= 0.513 #> kurtosis= 1.202 #> --------------------------------------------- #> #> Statistical tests #> Wilcoxon signed rank test : 0.652 #> Fisher variance test : 0.338 #> SW test of normality : 0.0757 . #> Global adjusted p-value : 0.227 #> --- #> Signif. codes: '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 #> ---------------------------------------------#> Performing simulations under the model. #> Plotting VPC #> Method used for VPC: binning by quantiles on X , dividing into the following intervals #> Interval Centered.On #> 1 (-1,3] 1.3 #> 2 (3,8] 7.4 #> 3 (8,14] 13.2 #> 4 (14,21] 20.5 #> 5 (21,37.7] 29.5 #> 6 (37.7,60] 50.4 #> 7 (60,90] 76.6 #> 8 (90,120] 109.0 #> 9 (120,180] 156.0f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:34 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:39 2020"#> Likelihoods computed by importance sampling#> AIC BIC #> 1 467.7644 465.0305 #> 2 469.4862 466.3617#>#>#># The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, # and compiled ODEs for FOMC, both are fast f_mmkin <- mmkin(list( "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), ds, quiet = TRUE) # These take about five seconds each on this system, as we use # analytical solutions written for saemix. When using the analytical # solutions written for mkin this took around four minutes f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:42 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:47 2020"f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])#> Running main SAEM algorithm #> [1] "Tue Nov 10 05:12:48 2020" #> .... #> Minimisation finished #> [1] "Tue Nov 10 05:12:57 2020"#> saemix version used for fitting: 3.1.9000 #> mkin version used for pre-fitting: 0.9.50.4 #> R version used for fitting: 4.0.3 #> Date of fit: Tue Nov 10 05:12:58 2020 #> Date of summary: Tue Nov 10 05:12:58 2020 #> #> Equations: #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * #> time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) #> * parent #> d_A1/dt = + f_parent_to_A1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) #> * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * #> exp(-k2 * time))) * parent - k_A1 * A1 #> #> Data: #> 170 observations of 2 variable(s) grouped in 5 datasets #> #> Model predictions using solution type analytical #> #> Fitted in 10.382 s using 300, 100 iterations #> #> Variance model: Constant variance #> #> Mean of starting values for individual parameters: #> parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 #> 93.8101519 -9.7647455 -0.9711148 -1.8799371 -4.2708142 #> g_qlogis #> 0.1356441 #> #> Fixed degradation parameter values: #> None #> #> Results: #> #> Likelihood computed by importance sampling #> AIC BIC logLik #> 841.3 836.2 -407.7 #> #> Optimised, transformed parameters with symmetric confidence intervals: #> est. lower upper #> parent_0 93.7514328 91.114 96.389 #> log_k_A1 -6.1262333 -8.432 -3.820 #> f_parent_qlogis -0.9739852 -1.372 -0.576 #> log_k1 -2.4818389 -3.747 -1.217 #> log_k2 -3.6138617 -5.294 -1.934 #> g_qlogis -0.0004614 -1.063 1.062 #> #> Correlation: #> prnt_0 lg__A1 f_prn_ log_k1 log_k2 #> log_k_A1 -0.013 #> f_parent_qlogis -0.025 0.050 #> log_k1 0.030 0.000 -0.005 #> log_k2 0.013 0.005 -0.003 0.037 #> g_qlogis -0.068 -0.016 0.011 -0.181 -0.181 #> #> Random effects: #> est. lower upper #> SD.parent_0 2.7857 0.7825 4.7889 #> SD.log_k_A1 2.1413 0.4425 3.8400 #> SD.f_parent_qlogis 0.4463 0.1609 0.7317 #> SD.log_k1 1.4097 0.5241 2.2954 #> SD.log_k2 1.8739 0.6979 3.0499 #> SD.g_qlogis 0.4559 -0.8150 1.7268 #> #> Variance model: #> est. lower upper #> a.1 1.883 1.666 2.1 #> #> Backtransformed parameters with asymmetric confidence intervals: #> est. lower upper #> parent_0 93.751433 9.111e+01 96.38921 #> k_A1 0.002185 2.177e-04 0.02193 #> f_parent_to_A1 0.274087 2.023e-01 0.35986 #> k1 0.083589 2.359e-02 0.29618 #> k2 0.026948 5.021e-03 0.14463 #> g 0.499885 2.567e-01 0.74312 #> #> Resulting formation fractions: #> ff #> parent_A1 0.2741 #> parent_sink 0.7259 #> #> Estimated disappearance times: #> DT50 DT90 DT50back DT50_k1 DT50_k2 #> parent 13.91 60.89 18.33 8.292 25.72 #> A1 317.26 1053.91 NA NA NA# Using a single core, the following takes about 6 minutes, using 10 cores # it is slower instead of faster #f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1) # }