Fit nonlinear mixed models with SAEM

mkinpredict( - x, - odeparms, - odeini, - outtimes = seq(0, 120, by = 0.1), - solution_type = "deSolve", - use_compiled = "auto", - method.ode = "lsoda", - atol = 1e-08, - rtol = 1e-10, - map_output = TRUE, - ... -) +

mkinpredict(x, odeparms, odeini, outtimes, ...)
 
 # S3 method for mkinmod
-mkinpredict(
-  x,
-  odeparms = c(k_parent_sink = 0.1),
-  odeini = c(parent = 100),
-  outtimes = seq(0, 120, by = 0.1),
-  solution_type = "deSolve",
-  use_compiled = "auto",
-  method.ode = "lsoda",
-  atol = 1e-08,
-  rtol = 1e-10,
-  map_output = TRUE,
-  ...
-)
+mkinpredict(
+  x,
+  odeparms = c(k_parent_sink = 0.1),
+  odeini = c(parent = 100),
+  outtimes = seq(0, 120, by = 0.1),
+  solution_type = "deSolve",
+  use_compiled = "auto",
+  method.ode = "lsoda",
+  atol = 1e-08,
+  rtol = 1e-10,
+  map_output = TRUE,
+  na_stop = TRUE,
+  ...
+)
 
 # S3 method for mkinfit
-mkinpredict(
-  x,
-  odeparms = x$bparms.ode,
-  odeini = x$bparms.state,
-  outtimes = seq(0, 120, by = 0.1),
-  solution_type = "deSolve",
-  use_compiled = "auto",
-  method.ode = "lsoda",
-  atol = 1e-08,
-  rtol = 1e-10,
-  map_output = TRUE,
-  ...
-)

+mkinpredict( + x, + odeparms = x$bparms.ode, + odeini = x$bparms.state, + outtimes = seq(0, 120, by = 0.1), + solution_type = "deSolve", + use_compiled = "auto", + method.ode = "lsoda", + atol = 1e-08, + rtol = 1e-10, + map_output = TRUE, + ... +)

Arguments

outtimes	A numeric vector specifying the time points for which model predictions should be generated.
...	Further arguments passed to the ode solver in case such a +solver is used.
solution_type
...	Further arguments passed to the ode solver in case such a -solver is used.	na_stop	Should it be an error if deSolve::ode returns NaN values

outtimes

A numeric vector specifying the time points for which model predictions should be generated.

...

Further arguments passed to the ode solver in case such a +solver is used.

solution_type

...

Further arguments passed to the ode solver in case such a -solver is used.

na_stop

Should it be an error if deSolve::ode returns NaN values

Examples

-SFO <- mkinmod(degradinol = mkinsub("SFO")) +SFO <- mkinmod(degradinol = mkinsub("SFO")) # Compare solution types -mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - solution_type = "analytical")

#> time degradinol +mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + solution_type = "analytical") +

#> time degradinol #> 0 0 100.0000000 #> 1 1 74.0818221 #> 2 2 54.8811636 @@ -294,8 +291,9 @@ solver is used.

#> 17 17 0.6096747 #> 18 18 0.4516581 #> 19 19 0.3345965 -#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - solution_type = "deSolve")

#> time degradinol +#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + solution_type = "deSolve") +

#> time degradinol #> 0 0 100.0000000 #> 1 1 74.0818221 #> 2 2 54.8811636 @@ -316,8 +314,9 @@ solver is used.

#> 17 17 0.6096747 #> 18 18 0.4516581 #> 19 19 0.3345965 -#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - solution_type = "deSolve", use_compiled = FALSE)

#> time degradinol +#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + solution_type = "deSolve", use_compiled = FALSE) +

#> time degradinol #> 0 0 100.0000000 #> 1 1 74.0818221 #> 2 2 54.8811636 @@ -338,8 +337,9 @@ solver is used.

#> 17 17 0.6096747 #> 18 18 0.4516581 #> 19 19 0.3345965 -#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - solution_type = "eigen")

#> time degradinol +#> 20 20 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + solution_type = "eigen") +

#> time degradinol #> 0 0 100.0000000 #> 1 1 74.0818221 #> 2 2 54.8811636 @@ -362,59 +362,68 @@ solver is used.

#> 19 19 0.3345965 #> 20 20 0.2478752

# Compare integration methods to analytical solution -mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - solution_type = "analytical")[21,]

#> time degradinol -#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - method = "lsoda")[21,]

#> time degradinol -#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - method = "ode45")[21,]

#> time degradinol -#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, - method = "rk4")[21,]

#> time degradinol +mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + solution_type = "analytical")[21,] +

#> time degradinol +#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + method = "lsoda")[21,] +

#> time degradinol +#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + method = "ode45")[21,] +

#> time degradinol +#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, + method = "rk4")[21,] +

#> time degradinol #> 20.0000000 0.2480043

# rk4 is not as precise here # The number of output times used to make a lot of difference until the # default for atol was adjusted -mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), - seq(0, 20, by = 0.1))[201,]

#> time degradinol -#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), - seq(0, 20, by = 0.01))[2001,]

#> time degradinol +mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), + seq(0, 20, by = 0.1))[201,] +

#> time degradinol +#> 20.0000000 0.2478752

mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), + seq(0, 20, by = 0.01))[2001,] +

#> time degradinol #> 20.0000000 0.2478752

# Comparison of the performance of solution types -SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), - m1 = list(type = "SFO"), use_of_ff = "max")

#> Successfully compiled differential equation model from auto-generated C code.

if(require(rbenchmark)) { - benchmark(replications = 10, order = "relative", columns = c("test", "relative", "elapsed"), - eigen = mkinpredict(SFO_SFO, - c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), - c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "eigen")[201,], - deSolve_compiled = mkinpredict(SFO_SFO, - c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), - c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "deSolve")[201,], - deSolve = mkinpredict(SFO_SFO, - c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), - c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "deSolve", use_compiled = FALSE)[201,], - analytical = mkinpredict(SFO_SFO, - c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), - c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "analytical", use_compiled = FALSE)[201,]) -}

#> test relative elapsed +SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), + m1 = list(type = "SFO"), use_of_ff = "max") +

#> Successfully compiled differential equation model from auto-generated C code.

if(require(rbenchmark)) { + benchmark(replications = 10, order = "relative", columns = c("test", "relative", "elapsed"), + eigen = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "eigen")[201,], + deSolve_compiled = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "deSolve")[201,], + deSolve = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "deSolve", use_compiled = FALSE)[201,], + analytical = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "analytical", use_compiled = FALSE)[201,]) +} +

#> Loading required package: rbenchmark

#> test relative elapsed #> 2 deSolve_compiled 1.0 0.005 -#> 4 analytical 1.8 0.009 #> 1 eigen 4.0 0.020 -#> 3 deSolve 44.6 0.223

+#> 4 analytical 4.0 0.020 +#> 3 deSolve 45.6 0.228

# \dontrun{ # Predict from a fitted model - f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) - f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "deSolve") - head(mkinpredict(f))

#> time parent m1 + f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) + f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "deSolve") + head(mkinpredict(f)) +

#> time parent m1 #> 0 0.0 82.49216 0.000000 -#> 0.1 0.1 80.00562 1.236198 -#> 0.2 0.2 77.59404 2.422818 -#> 0.3 0.3 75.25514 3.561476 -#> 0.4 0.4 72.98675 4.653740 -#> 0.5 0.5 70.78673 5.701130

# } +#> 0.1 0.1 80.00562 1.236394 +#> 0.2 0.2 77.59404 2.423201 +#> 0.3 0.3 75.25514 3.562040 +#> 0.4 0.4 72.98675 4.654478 +#> 0.5 0.5 70.78673 5.702033

# }

Site built with pkgdown 1.5.1.

Site built with pkgdown 1.6.1.

Fit nonlinear mixed models with SAEM — saem • mkin

+ + + + +

+ +

+ + +

This function uses saemix::saemix() as a backend for fitting nonlinear mixed +effects models created from mmkin row objects using the stochastic approximation +to the 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,
+  ...
+)
+
+saemix_model(object, cores = 1, verbose = FALSE, ...)
+
+saemix_data(object, verbose = FALSE, ...)

+ +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + + + + + +

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 +`parallel::mclapply()`. Using more than 1 core is experimental and may +lead to uncontrolled forking, apparently depending on the BLAS version +used.
verbose	Should we print information about created objects?
suppressPlot	Should we suppress any plotting that is done +by the saemix function?

+ +

Value

+ +

An saemix::SaemixObject.

An saemix::SaemixModel object.

An saemix::SaemixData object.

Details

+ +

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.

+ +

Examples

# \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] "Sat Nov  7 13:14:50 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:14:52 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] "Sat Nov  7 13:14:53 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:14:55 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+
#> Running main SAEM algorithm
+#> [1] "Sat Nov  7 13:14:55 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:14:57 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+
#> Running main SAEM algorithm
+#> [1] "Sat Nov  7 13:14:57 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:15:00 2020"

+# The returned saem.mmkin object contains an SaemixObject, we can use
+# functions from saemix
+library(saemix)
+
#> Package saemix, version 3.1.9000
+#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
+
#> Likelihoods computed by importance sampling 
#>        AIC      BIC
+#> 1 624.2428 622.2900
+#> 2 467.7644 465.0305
+#> 3 491.3541 487.8391

+f_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] "Sat Nov  7 13:15:02 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:15:07 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+
#> Likelihoods computed by importance sampling 
#>        AIC      BIC
+#> 1 467.7644 465.0305
+#> 2 469.4862 466.3617

+dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+  A1 = mkinsub("SFO"))
+
#> Successfully compiled differential equation model from auto-generated C code.
f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
+# This takes about 4 minutes on my system
+f_saem <- saem(f_mmkin)
+
#> Running main SAEM algorithm
+#> [1] "Sat Nov  7 13:15:08 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:19:07 2020"

+f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
+# Using a single core, the following takes about 6 minutes, using 10 cores
+# it is slower instead of faster
+f_saem_des <- saem(f_mmkin_des, cores = 1)
+
#> Running main SAEM algorithm
+#> [1] "Sat Nov  7 13:19:26 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Sat Nov  7 13:27:33 2020"
compare.saemix(list(f_saemix$so, f_saemix_des$so))
+
#> Error in compare.saemix(list(f_saemix$so, f_saemix_des$so)): object 'f_saemix' not found
# }
+

+ +

+ + + +

# S3 method for mkinfit -summary(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) +summary(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) # S3 method for summary.mkinfit -print(x, digits = max(3, getOption("digits") - 3), ...)

Examples

- summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE))

#> mkin version used for fitting: 0.9.50.3 -#> R version used for fitting: 4.0.2 -#> Date of fit: Thu Oct 8 09:13:59 2020 -#> Date of summary: Thu Oct 8 09:13:59 2020 + summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE)) +

#> mkin version used for fitting: 0.9.50.4 +#> R version used for fitting: 4.0.3 +#> Date of fit: Sat Nov 7 13:28:04 2020 +#> Date of summary: Sat Nov 7 13:28:04 2020 #> #> Equations: #> d_parent/dt = - k_parent * parent #> #> Model predictions using solution type analytical #> -#> Fitted using 131 model solutions performed in 0.027 s +#> Fitted using 131 model solutions performed in 0.028 s #> #> Error model: Constant variance #> @@ -322,7 +326,7 @@ EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

Site built with pkgdown 1.5.1.

Site built with pkgdown 1.6.1.

diff --git a/docs/dev/reference/summary.nlme.mmkin.html b/docs/dev/reference/summary.nlme.mmkin.html index 38f53122..d35e8713 100644 --- a/docs/dev/reference/summary.nlme.mmkin.html +++ b/docs/dev/reference/summary.nlme.mmkin.html @@ -156,7 +156,7 @@ endpoints such as formation fractions and DT50 values. Optionally

# S3 method for nlme.mmkin
-summary(
+summary(
   object,
   data = FALSE,
   verbose = FALSE,
@@ -260,13 +260,13 @@ José Pinheiro and Douglas Bates for the components inherited from nlme
 # Evaluate using mmkin and nlme
 library(nlme)
 f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
-#> Warning: Shapiro-Wilk test for standardized residuals: p =  0.00371
f_nlme <- nlme(f_mmkin)
-
#> Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
summary(f_nlme, data = TRUE)
+f_nlme <- nlme(f_mmkin)
+
#> Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
summary(f_nlme, data = TRUE)
 
#> nlme version used for fitting:      3.1.150.1 
 #> mkin version used for pre-fitting:  0.9.50.4 
 #> R version used for fitting:         4.0.3 
-#> Date of fit:     Tue Oct 27 15:25:23 2020 
-#> Date of summary: Tue Oct 27 15:25:23 2020 
+#> Date of fit:     Sat Nov  7 13:28:06 2020 
+#> Date of summary: Sat Nov  7 13:28:06 2020 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -276,7 +276,7 @@ José Pinheiro and Douglas Bates for the components inherited from nlme
 #> 
 #> Model predictions using solution type analytical 
 #> 
-#> Fitted in 0.581 s using 5 iterations
+#> Fitted in 0.588 s using 5 iterations
 #> 
 #> Variance model: Two-component variance function 
 #> 
-- 
cgit v1.2.1

Produce predictions from a kinetic model using specific parameters

Arguments

Value

Author

Examples

Fit nonlinear mixed models with SAEM

Arguments

Value

Details

Examples

Summary method for class "mkinfit"

Arguments

Author

Examples