From c6a1733974334b4e97a27170c60e481dc9e9f35d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 11 Dec 2020 15:39:43 +0100 Subject: Update static docs --- docs/dev/reference/Rplot001.png | Bin 14180 -> 19395 bytes docs/dev/reference/Rplot002.png | Bin 13753 -> 16843 bytes docs/dev/reference/Rplot005.png | Bin 56740 -> 57095 bytes docs/dev/reference/index.html | 6 + docs/dev/reference/mixed-1.png | Bin 0 -> 206373 bytes docs/dev/reference/mixed.html | 279 ++++++++++++++++++++++++++++++ docs/dev/reference/mkinmod.html | 8 +- docs/dev/reference/nlme.mmkin-1.png | Bin 122068 -> 119655 bytes docs/dev/reference/nlme.mmkin-2.png | Bin 160690 -> 159253 bytes docs/dev/reference/nlme.mmkin.html | 106 +++++------- docs/dev/reference/plot.mixed.mmkin-2.png | Bin 163684 -> 164488 bytes docs/dev/reference/plot.mixed.mmkin-3.png | Bin 162677 -> 163488 bytes docs/dev/reference/plot.mixed.mmkin.html | 8 +- docs/dev/reference/saem-5.png | Bin 163636 -> 164443 bytes docs/dev/reference/saem.html | 127 +++++++------- 15 files changed, 399 insertions(+), 135 deletions(-) create mode 100644 docs/dev/reference/mixed-1.png create mode 100644 docs/dev/reference/mixed.html (limited to 'docs/dev/reference') diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png index ae5db971..bca41e2c 100644 Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ diff --git a/docs/dev/reference/Rplot002.png b/docs/dev/reference/Rplot002.png index 51e51cf4..9b97a634 100644 Binary files a/docs/dev/reference/Rplot002.png and b/docs/dev/reference/Rplot002.png differ diff --git a/docs/dev/reference/Rplot005.png b/docs/dev/reference/Rplot005.png index 91c7777b..5e675828 100644 Binary files a/docs/dev/reference/Rplot005.png and b/docs/dev/reference/Rplot005.png differ diff --git a/docs/dev/reference/index.html b/docs/dev/reference/index.html index cb37f9a6..5c964cf9 100644 --- a/docs/dev/reference/index.html +++ b/docs/dev/reference/index.html @@ -358,6 +358,12 @@ of an mmkin object

get_deg_func()

Retrieve a degradation function from the mmkin namespace

+ + + +

mixed() print(<mixed.mmkin>)

+ +

Create a mixed effects model from an mmkin row object

diff --git a/docs/dev/reference/mixed-1.png b/docs/dev/reference/mixed-1.png new file mode 100644 index 00000000..3400c4aa Binary files /dev/null and b/docs/dev/reference/mixed-1.png differ diff --git a/docs/dev/reference/mixed.html b/docs/dev/reference/mixed.html new file mode 100644 index 00000000..18a67af8 --- /dev/null +++ b/docs/dev/reference/mixed.html @@ -0,0 +1,279 @@ + + + + + + + + +Create a mixed effects model from an mmkin row object — mixed • mkin + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Create a mixed effects model from an mmkin row object

+ Source: R/mixed.mmkin.R + +
+ +
+

Create a mixed effects model from an mmkin row object

+
+ + 
mixed(object, ...)
+
+# S3 method for mmkin
+mixed(object, method = c("none"), ...)
+
+# S3 method for mixed.mmkin
+print(x, digits = max(3, getOption("digits") - 3), ...)
+ +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + +
object

An mmkin row object

...

Currently not used

method

The method to be used

x

A mixed.mmkin object to print

digits

Number of digits to use for printing.

+ + +

Examples

+ 
sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+n_biphasic <- 8
+err_1 = list(const = 1, prop = 0.07)
+
+DFOP_SFO <- mkinmod(
+  parent = mkinsub("DFOP", "m1"),
+  m1 = mkinsub("SFO"),
+  quiet = TRUE)
+
+set.seed(123456)
+log_sd <- 0.3
+syn_biphasic_parms <- as.matrix(data.frame(
+  k1 = rlnorm(n_biphasic, log(0.05), log_sd),
+  k2 = rlnorm(n_biphasic, log(0.01), log_sd),
+  g = plogis(rnorm(n_biphasic, 0, log_sd)),
+  f_parent_to_m1 = plogis(rnorm(n_biphasic, 0, log_sd)),
+  k_m1 = rlnorm(n_biphasic, log(0.002), log_sd)))
+
+ds_biphasic_mean <- lapply(1:n_biphasic,
+  function(i) {
+    mkinpredict(DFOP_SFO, syn_biphasic_parms[i, ],
+      c(parent = 100, m1 = 0), sampling_times)
+  }
+)
+
+set.seed(123456L)
+ds_biphasic <- lapply(ds_biphasic_mean, function(ds) {
+  add_err(ds,
+    sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2),
+    n = 1, secondary = "m1")[[1]]
+})
+
+# \dontrun{
+f_mmkin <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, error_model = "tc", quiet = TRUE)
+
+f_mixed <- mixed(f_mmkin)
+print(f_mixed)
+
#> Kinetic model fitted by nonlinear regression to each dataset
+#> Structural model:
+#> 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_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g)
+#>            * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) *
+#>            exp(-k2 * time))) * parent - k_m1 * m1
+#> 
+#> Data:
+#> 271 observations of 2 variable(s) grouped in 8 datasets
+#> 
+#> <mmkin> object
+#> Status of individual fits:
+#> 
+#>           dataset
+#> model      1  2  3  4  5  6 7  8 
+#>   DFOP-SFO OK OK OK OK OK C OK OK
+#> 
+#> OK: No warnings
+#> C: Optimisation did not converge:
+#> iteration limit reached without convergence (10)
+#> 
+#> Mean fitted parameters:
+#>        parent_0        log_k_m1 f_parent_qlogis          log_k1          log_k2 
+#>      100.606304       -8.759216       -0.002001       -3.350539       -3.989549 
+#>        g_qlogis 
+#>       -0.090353 
plot(f_mixed)
+
# }
+
+
+ +
+ + +
+
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 1.6.1.

+
+ +
+
+ + + + + + + + diff --git a/docs/dev/reference/mkinmod.html b/docs/dev/reference/mkinmod.html index 9e37e664..9a7dac6f 100644 --- a/docs/dev/reference/mkinmod.html +++ b/docs/dev/reference/mkinmod.html @@ -201,8 +201,8 @@ In print.mkinmod, this argument is currently not used.

Specification of the use of formation fractions in the model equations and, if applicable, the coefficient matrix. If "max", formation fractions are always used (default). If "min", a minimum use of -formation fractions is made, i.e. each pathway to a metabolite has its -own rate constant.

+formation fractions is made, i.e. each first-order pathway to a metabolite +has its own rate constant.

name @@ -348,7 +348,7 @@ Evaluating and Calculating Degradation Kinetics in Environmental Media

parent = mkinsub("SFO", "m1", full_name = "Test compound"), m1 = mkinsub("SFO", full_name = "Metabolite M1"), name = "SFO_SFO", dll_dir = DLL_dir, unload = TRUE, overwrite = TRUE) -
#> Copied DLL from /tmp/RtmpX9Xf3y/file6d622f3b7e08.so to /home/jranke/.local/share/mkin/SFO_SFO.so
# Now we can save the model and restore it in a new session +
#> Copied DLL from /tmp/Rtmpy4eiQb/file554e573761a7.so to /home/jranke/.local/share/mkin/SFO_SFO.so
# Now we can save the model and restore it in a new session saveRDS(SFO_SFO.2, file = "~/SFO_SFO.rds") # Terminate the R session here if you would like to check, and then do library(mkin) @@ -397,7 +397,7 @@ Evaluating and Calculating Degradation Kinetics in Environmental Media

#> }) #> return(predicted) #> } -#> <environment: 0x55555c5b8d68>
+#> <environment: 0x55555cac8d00>
# If we have several parallel metabolites # (compare tests/testthat/test_synthetic_data_for_UBA_2014.R) m_synth_DFOP_par <- mkinmod( diff --git a/docs/dev/reference/nlme.mmkin-1.png b/docs/dev/reference/nlme.mmkin-1.png index d760e8f9..25bebeca 100644 Binary files a/docs/dev/reference/nlme.mmkin-1.png and b/docs/dev/reference/nlme.mmkin-1.png differ diff --git a/docs/dev/reference/nlme.mmkin-2.png b/docs/dev/reference/nlme.mmkin-2.png index c7085b81..c314c149 100644 Binary files a/docs/dev/reference/nlme.mmkin-2.png and b/docs/dev/reference/nlme.mmkin-2.png differ diff --git a/docs/dev/reference/nlme.mmkin.html b/docs/dev/reference/nlme.mmkin.html index 6d9f2007..a4d7070a 100644 --- a/docs/dev/reference/nlme.mmkin.html +++ b/docs/dev/reference/nlme.mmkin.html @@ -154,11 +154,12 @@ have been obtained by fitting the same model to a list of datasets.

# S3 method for mmkin
 nlme(
   model,
-  data = sys.frame(sys.parent()),
-  fixed,
-  random = fixed,
+  data = "auto",
+  fixed = lapply(as.list(names(mean_degparms(model))), function(el) eval(parse(text =
+    paste(el, 1, sep = "~")))),
+  random = pdDiag(fixed),
   groups,
-  start,
+  start = mean_degparms(model, random = TRUE),
   correlation = NULL,
   weights = NULL,
   subset,
@@ -276,14 +277,14 @@ methods that will automatically work on 'nlme.mmkin' objects, such as
 f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1)
 library(nlme)
 f_nlme_sfo <- nlme(f["SFO", ])
-
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
+ # \dontrun{ f_nlme_dfop <- nlme(f["DFOP", ]) anova(f_nlme_sfo, f_nlme_dfop)
#> Model df AIC BIC logLik Test L.Ratio p-value -#> f_nlme_sfo 1 6 622.0677 637.0666 -305.0338 -#> f_nlme_dfop 2 15 487.0134 524.5105 -228.5067 1 vs 2 153.0543 <.0001
print(f_nlme_dfop) +#> f_nlme_sfo 1 5 625.0539 637.5529 -307.5269 +#> f_nlme_dfop 2 9 495.1270 517.6253 -238.5635 1 vs 2 137.9268 <.0001
print(f_nlme_dfop)
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood #> #> Structural model: @@ -294,28 +295,24 @@ methods that will automatically work on 'nlme.mmkin' objects, such as #> Data: #> 90 observations of 1 variable(s) grouped in 5 datasets #> -#> Log-likelihood: -228.5067 +#> Log-likelihood: -238.6 #> #> Fixed effects: #> list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) #> parent_0 log_k1 log_k2 g_qlogis -#> 94.18273 -1.82135 -4.16872 0.08949 +#> 94.1702 -1.8002 -4.1474 0.0324 #> #> Random effects: #> Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) #> Level: ds -#> Structure: General positive-definite, Log-Cholesky parametrization -#> StdDev Corr -#> parent_0 2.4656397 prnt_0 log_k1 log_k2 -#> log_k1 0.7950788 0.240 -#> log_k2 1.2605419 0.150 0.984 -#> g_qlogis 0.5013272 -0.075 0.843 0.834 -#> Residual 2.3308100 +#> Structure: Diagonal +#> parent_0 log_k1 log_k2 g_qlogis Residual +#> StdDev: 2.488 0.8447 1.33 0.4652 2.321 #>
plot(f_nlme_dfop)
endpoints(f_nlme_dfop)
#> $distimes #> DT50 DT90 DT50back DT50_k1 DT50_k2 -#> parent 10.57119 101.0652 30.42366 4.283776 44.80015 +#> parent 10.79857 100.7937 30.34192 4.193937 43.85442 #>
ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], function(x) x$data[c("name", "time", "value")]) @@ -332,7 +329,7 @@ methods that will automatically work on 'nlme.mmkin' objects, such as ds_2, quiet = TRUE) f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ]) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
plot(f_nlme_sfo_sfo) + plot(f_nlme_sfo_sfo)
# With formation fractions this does not coverge with defaults # f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) @@ -343,32 +340,22 @@ methods that will automatically work on 'nlme.mmkin' objects, such as # parameters (with pdDiag, increasing the tolerance and pnlsMaxIter was # necessary) f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ]) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 4, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
+
#> Error in nlme.formula(model = value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_A1, f_parent_qlogis, log_k1, log_k2, g_qlogis), data = structure(list(ds = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L ), .Label = c("1", "2", "3", "4", "5"), class = c("ordered", "factor")), name = c("parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1"), time = c(0, 0, 3, 3, 6, 6, 10, 10, 20, 20, 34, 34, 55, 55, 90, 90, 112, 112, 132, 132, 3, 3, 6, 6, 10, 10, 20, 20, 34, 34, 55, 55, 90, 90, 112, 112, 132, 132, 0, 0, 3, 3, 7, 7, 14, 14, 30, 30, 60, 60, 90, 90, 120, 120, 180, 180, 3, 3, 7, 7, 14, 14, 30, 30, 60, 60, 90, 90, 120, 120, 180, 180, 0, 0, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 0, 0, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 91, 91, 120, 120, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 91, 91, 120, 120, 0, 0, 8, 8, 14, 14, 21, 21, 41, 41, 63, 63, 91, 91, 120, 120, 8, 8, 14, 14, 21, 21, 41, 41, 63, 63, 91, 91, 120, 120), value = c(97.2, 96.4, 71.1, 69.2, 58.1, 56.6, 44.4, 43.4, 33.3, 29.2, 17.6, 18, 10.5, 9.3, 4.5, 4.7, 3, 3.4, 2.3, 2.7, 4.3, 4.6, 7, 7.2, 8.2, 8, 11, 13.7, 11.5, 12.7, 14.9, 14.5, 12.1, 12.3, 9.9, 10.2, 8.8, 7.8, 93.6, 92.3, 87, 82.2, 74, 73.9, 64.2, 69.5, 54, 54.6, 41.1, 38.4, 32.5, 35.5, 28.1, 29, 26.5, 27.6, 3.9, 3.1, 6.9, 6.6, 10.4, 8.3, 14.4, 13.7, 22.1, 22.3, 27.5, 25.4, 28, 26.6, 25.8, 25.3, 91.9, 90.8, 64.9, 66.2, 43.5, 44.1, 18.3, 18.1, 10.2, 10.8, 4.9, 3.3, 1.6, 1.5, 1.1, 0.9, 9.6, 7.7, 15, 15.1, 21.2, 21.1, 19.7, 18.9, 17.5, 15.9, 9.5, 9.8, 6.2, 6.1, 99.8, 98.3, 77.1, 77.2, 59, 58.1, 27.4, 29.2, 19.1, 29.6, 10.1, 18.2, 4.5, 9.1, 2.3, 2.9, 2, 1.8, 2, 2.2, 4.2, 3.9, 7.4, 7.9, 14.5, 13.7, 14.2, 12.2, 13.7, 13.2, 13.6, 15.4, 10.4, 11.6, 10, 9.5, 9.1, 9, 96.1, 94.3, 73.9, 73.9, 69.4, 73.1, 65.6, 65.3, 55.9, 54.4, 47, 49.3, 44.7, 46.7, 42.1, 41.3, 3.3, 3.4, 3.9, 2.9, 6.4, 7.2, 9.1, 8.5, 11.7, 12, 13.3, 13.2, 14.3, 12.1)), row.names = c(NA, -170L), class = c("nfnGroupedData", "nfGroupedData", "groupedData", "data.frame"), formula = value ~ time | ds, FUN = function (x) max(x, na.rm = TRUE), order.groups = FALSE), start = list( fixed = c(parent_0 = 93.8101519326534, log_k_A1 = -9.76474551635931, f_parent_qlogis = -0.971114801595408, log_k1 = -1.87993711571859, log_k2 = -4.27081421366622, g_qlogis = 0.135644115277507 ), random = list(ds = structure(c(2.56569977430371, -3.49441920289139, -3.32614443321494, 4.35347873814922, -0.0986148763466161, 4.65850590018027, 1.8618544764481, 6.12693257601545, 4.91792724701579, -17.5652201996596, -0.466203822618637, 0.746660653597927, 0.282193987271096, -0.42053488943072, -0.142115928819667, 0.369240076779088, -1.38985563501659, 1.02592753494098, 0.73090914081534, -0.736221117518819, 0.768170629350299, -1.89347658079869, 1.72168783460352, 0.844607177798114, -1.44098906095325, -0.377731855445672, 0.168180098477565, 0.469683412912104, 0.500717664434525, -0.760849320378522), .Dim = 5:6, .Dimnames = list(c("1", "2", "3", "4", "5"), c("parent_0", "log_k_A1", "f_parent_qlogis", "log_k1", "log_k2", "g_qlogis"))))), fixed = list(parent_0 ~ 1, log_k_A1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1), random = structure(numeric(0), class = c("pdDiag", "pdMat"), formula = structure(list(parent_0 ~ 1, log_k_A1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1), class = "listForm"), Dimnames = list(NULL, NULL))): maximum number of iterations (maxIter = 50) reached without convergence
#> Timing stopped at: 49.95 16.5 44.08
plot(f_nlme_dfop_sfo) -
+
#> Error in plot(f_nlme_dfop_sfo): object 'f_nlme_dfop_sfo' not found
anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) -
#> Model df AIC BIC logLik Test L.Ratio p-value -#> f_nlme_dfop_sfo 1 28 811.7199 899.5222 -377.8599 -#> f_nlme_sfo_sfo 2 15 1075.1934 1122.2304 -522.5967 1 vs 2 289.4736 <.0001
+
#> Error in anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo): object 'f_nlme_dfop_sfo' not found
endpoints(f_nlme_sfo_sfo)
#> $ff #> parent_sink parent_A1 A1_sink -#> 0.6512742 0.3487258 1.0000000 +#> 0.5912432 0.4087568 1.0000000 #> #> $distimes -#> DT50 DT90 -#> parent 18.03144 59.89916 -#> A1 102.72949 341.25997 +#> DT50 DT90 +#> parent 19.13518 63.5657 +#> A1 66.02155 219.3189 #>
endpoints(f_nlme_dfop_sfo) -
#> $ff -#> parent_A1 parent_sink -#> 0.2762167 0.7237833 -#> -#> $distimes -#> DT50 DT90 DT50back DT50_k1 DT50_k2 -#> parent 11.15024 133.9652 40.32755 4.688015 62.16017 -#> A1 235.83191 783.4167 NA NA NA -#>
+
#> Error in endpoints(f_nlme_dfop_sfo): object 'f_nlme_dfop_sfo' not found
if (length(findFunction("varConstProp")) > 0) { # tc error model for nlme available # Attempts to fit metabolite kinetics with the tc error model are possible, # but need tweeking of control values and sometimes do not converge @@ -379,7 +366,7 @@ methods that will automatically work on 'nlme.mmkin' objects, such as AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc) print(f_nlme_dfop_tc) } -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 14, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 4, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 5, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 6, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 7, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 8, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 9, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 10, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 11, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 12, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 14, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 15, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 16, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 17, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 18, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood +
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood #> #> Structural model: #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -389,35 +376,31 @@ methods that will automatically work on 'nlme.mmkin' objects, such as #> Data: #> 90 observations of 1 variable(s) grouped in 5 datasets #> -#> Log-likelihood: -228.3575 +#> Log-likelihood: -238.4 #> #> Fixed effects: #> list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) #> parent_0 log_k1 log_k2 g_qlogis -#> 93.6695 -1.9187 -4.4253 0.2215 +#> 94.04775 -1.82340 -4.16715 0.05685 #> #> Random effects: #> Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) #> Level: ds -#> Structure: General positive-definite, Log-Cholesky parametrization -#> StdDev Corr -#> parent_0 2.8574651 prnt_0 log_k1 log_k2 -#> log_k1 0.9689083 0.506 -#> log_k2 1.5798002 0.446 0.997 -#> g_qlogis 0.5761569 -0.457 0.247 0.263 -#> Residual 1.0000000 +#> Structure: Diagonal +#> parent_0 log_k1 log_k2 g_qlogis Residual +#> StdDev: 2.474 0.85 1.337 0.4659 1 #> #> Variance function: #> Structure: Constant plus proportion of variance covariate #> Formula: ~fitted(.) #> Parameter estimates: -#> const prop -#> 2.0376990 0.0221686
+#> const prop +#> 2.23224114 0.01262341
f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo, "DFOP-SFO" = m_dfop_sfo), ds_2, quiet = TRUE, error_model = "obs") f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
print(f_nlme_sfo_sfo_obs) + print(f_nlme_sfo_sfo_obs)
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood #> #> Structural model: @@ -427,44 +410,37 @@ methods that will automatically work on 'nlme.mmkin' objects, such as #> Data: #> 170 observations of 2 variable(s) grouped in 5 datasets #> -#> Log-likelihood: -462.2203 +#> Log-likelihood: -473 #> #> Fixed effects: #> list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1, log_k_A1_sink ~ 1) #> parent_0 log_k_parent_sink log_k_parent_A1 log_k_A1_sink -#> 88.682 -3.664 -4.164 -4.665 +#> 87.976 -3.670 -4.164 -4.645 #> #> Random effects: #> Formula: list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1, log_k_A1_sink ~ 1) #> Level: ds -#> Structure: General positive-definite, Log-Cholesky parametrization -#> StdDev Corr -#> parent_0 4.9153305 prnt_0 lg_k__ l___A1 -#> log_k_parent_sink 1.8158570 0.956 -#> log_k_parent_A1 1.0514548 0.821 0.907 -#> log_k_A1_sink 0.4924122 0.035 0.315 0.533 -#> Residual 6.3987599 +#> Structure: Diagonal +#> parent_0 log_k_parent_sink log_k_parent_A1 log_k_A1_sink Residual +#> StdDev: 3.992 1.777 1.055 0.4821 6.483 #> #> Variance function: #> Structure: Different standard deviations per stratum #> Formula: ~1 | name #> Parameter estimates: #> parent A1 -#> 1.0000000 0.2040647
f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ]) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 4, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
+#> 1.0000000 0.2050003
f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ]) +
#> Error in nlme.formula(model = value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_A1, f_parent_qlogis, log_k1, log_k2, g_qlogis), data = structure(list(ds = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L ), .Label = c("1", "2", "3", "4", "5"), class = c("ordered", "factor")), name = c("parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "parent", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1"), time = c(0, 0, 3, 3, 6, 6, 10, 10, 20, 20, 34, 34, 55, 55, 90, 90, 112, 112, 132, 132, 3, 3, 6, 6, 10, 10, 20, 20, 34, 34, 55, 55, 90, 90, 112, 112, 132, 132, 0, 0, 3, 3, 7, 7, 14, 14, 30, 30, 60, 60, 90, 90, 120, 120, 180, 180, 3, 3, 7, 7, 14, 14, 30, 30, 60, 60, 90, 90, 120, 120, 180, 180, 0, 0, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 0, 0, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 91, 91, 120, 120, 1, 1, 3, 3, 8, 8, 14, 14, 27, 27, 48, 48, 70, 70, 91, 91, 120, 120, 0, 0, 8, 8, 14, 14, 21, 21, 41, 41, 63, 63, 91, 91, 120, 120, 8, 8, 14, 14, 21, 21, 41, 41, 63, 63, 91, 91, 120, 120), value = c(97.2, 96.4, 71.1, 69.2, 58.1, 56.6, 44.4, 43.4, 33.3, 29.2, 17.6, 18, 10.5, 9.3, 4.5, 4.7, 3, 3.4, 2.3, 2.7, 4.3, 4.6, 7, 7.2, 8.2, 8, 11, 13.7, 11.5, 12.7, 14.9, 14.5, 12.1, 12.3, 9.9, 10.2, 8.8, 7.8, 93.6, 92.3, 87, 82.2, 74, 73.9, 64.2, 69.5, 54, 54.6, 41.1, 38.4, 32.5, 35.5, 28.1, 29, 26.5, 27.6, 3.9, 3.1, 6.9, 6.6, 10.4, 8.3, 14.4, 13.7, 22.1, 22.3, 27.5, 25.4, 28, 26.6, 25.8, 25.3, 91.9, 90.8, 64.9, 66.2, 43.5, 44.1, 18.3, 18.1, 10.2, 10.8, 4.9, 3.3, 1.6, 1.5, 1.1, 0.9, 9.6, 7.7, 15, 15.1, 21.2, 21.1, 19.7, 18.9, 17.5, 15.9, 9.5, 9.8, 6.2, 6.1, 99.8, 98.3, 77.1, 77.2, 59, 58.1, 27.4, 29.2, 19.1, 29.6, 10.1, 18.2, 4.5, 9.1, 2.3, 2.9, 2, 1.8, 2, 2.2, 4.2, 3.9, 7.4, 7.9, 14.5, 13.7, 14.2, 12.2, 13.7, 13.2, 13.6, 15.4, 10.4, 11.6, 10, 9.5, 9.1, 9, 96.1, 94.3, 73.9, 73.9, 69.4, 73.1, 65.6, 65.3, 55.9, 54.4, 47, 49.3, 44.7, 46.7, 42.1, 41.3, 3.3, 3.4, 3.9, 2.9, 6.4, 7.2, 9.1, 8.5, 11.7, 12, 13.3, 13.2, 14.3, 12.1)), row.names = c(NA, -170L), class = c("nfnGroupedData", "nfGroupedData", "groupedData", "data.frame"), formula = value ~ time | ds, FUN = function (x) max(x, na.rm = TRUE), order.groups = FALSE), start = list( fixed = c(parent_0 = 93.4272167134207, log_k_A1 = -9.71590717106959, f_parent_qlogis = -0.953712099744438, log_k1 = -1.95256957646888, log_k2 = -4.42919226610318, g_qlogis = 0.193023137298073 ), random = list(ds = structure(c(2.85557330683041, -3.87630303729395, -2.78062140212751, 4.82042042600536, -1.01906929341432, 4.613992019697, 2.05871276943309, 6.0766404049189, 4.86471337131288, -17.6140585653619, -0.480721175257541, 0.773079218835614, 0.260464433006093, -0.440615012802434, -0.112207463781733, 0.445812953745225, -1.49588630006094, 1.13602040717272, 0.801850880762046, -0.887797941619048, 0.936480292463262, -2.43093808171905, 1.91256225793793, 0.984827519864443, -1.40293198854659, -0.455176326336681, 0.376355651864385, 0.343919720700401, 0.46329187713133, -0.728390923359434 ), .Dim = 5:6, .Dimnames = list(c("1", "2", "3", "4", "5"), c("parent_0", "log_k_A1", "f_parent_qlogis", "log_k1", "log_k2", "g_qlogis"))))), fixed = list(parent_0 ~ 1, log_k_A1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1), random = structure(numeric(0), class = c("pdDiag", "pdMat"), formula = structure(list(parent_0 ~ 1, log_k_A1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1), class = "listForm"), Dimnames = list(NULL, NULL)), weights = structure(numeric(0), formula = ~1 | name, class = c("varIdent", "varFunc"))): maximum number of iterations (maxIter = 50) reached without convergence
#> Timing stopped at: 59.38 16.5 53.5
f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, "DFOP-SFO" = m_dfop_sfo), ds_2, quiet = TRUE, error_model = "tc") # f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ]) # stops with error message f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 4, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 6, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 7, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 8, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 9, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 11, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 12, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 15, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
#> Warning: Iteration 25, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
# We get warnings about false convergence in the LME step in several iterations +
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Warning: longer object length is not a multiple of shorter object length
#> Error in X[, fmap[[nm]]] <- gradnm: number of items to replace is not a multiple of replacement length
#> Timing stopped at: 6.363 2.688 5.469
# We get warnings about false convergence in the LME step in several iterations # but as the last such warning occurs in iteration 25 and we have 28 iterations # we can ignore these anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs, f_nlme_dfop_sfo_tc) -
#> Model df AIC BIC logLik Test L.Ratio p-value -#> f_nlme_dfop_sfo 1 28 811.7199 899.5222 -377.8599 -#> f_nlme_dfop_sfo_obs 2 29 784.1304 875.0685 -363.0652 1 vs 2 29.5895 <.0001 -#> f_nlme_dfop_sfo_tc 3 29 791.9981 882.9362 -366.9990
+
#> Error in anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs, f_nlme_dfop_sfo_tc): object 'f_nlme_dfop_sfo' not found
# }
diff --git a/docs/dev/reference/plot.mixed.mmkin-2.png b/docs/dev/reference/plot.mixed.mmkin-2.png index 0821dcca..c0d67204 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-2.png and b/docs/dev/reference/plot.mixed.mmkin-2.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin-3.png b/docs/dev/reference/plot.mixed.mmkin-3.png index d58b9c69..3055c0c9 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-3.png and b/docs/dev/reference/plot.mixed.mmkin-3.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin.html b/docs/dev/reference/plot.mixed.mmkin.html index 492e289c..90ba9184 100644 --- a/docs/dev/reference/plot.mixed.mmkin.html +++ b/docs/dev/reference/plot.mixed.mmkin.html @@ -174,7 +174,7 @@ x -

An object of class saem.mmkin or nlme.mmkin

+

An object of class mixed.mmkin, saem.mmkin or nlme.mmkin

i @@ -273,14 +273,14 @@ corresponding model prediction lines for the different datasets.

# For this fit we need to increase pnlsMaxiter, and we increase the # tolerance in order to speed up the fit for this example evaluation f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) -
#> Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
#> Warning: Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!
plot(f_nlme) +plot(f_nlme)
f_saem <- saem(f)
#> Running main SAEM algorithm -#> [1] "Mon Nov 30 15:52:45 2020" +#> [1] "Fri Dec 11 15:37:36 2020" #> .... #> Minimisation finished -#> [1] "Mon Nov 30 15:52:54 2020"
plot(f_saem) +#> [1] "Fri Dec 11 15:37:45 2020"
plot(f_saem)
# }
diff --git a/docs/dev/reference/saem-5.png b/docs/dev/reference/saem-5.png index a12f8383..6e6e0f91 100644 Binary files a/docs/dev/reference/saem-5.png and b/docs/dev/reference/saem-5.png differ diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html index e4ebf5d9..73699e19 100644 --- a/docs/dev/reference/saem.html +++ b/docs/dev/reference/saem.html @@ -151,11 +151,13 @@ effects models created from mmkin row objects using the Expectation Maximisation algorithm (SAEM).

- 
saem(object, control, ...)
+    
saem(object, ...)
 
 # S3 method for mmkin
 saem(
   object,
+  transformations = c("mkin", "saemix"),
+  solution_type = "auto",
   control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs =
     FALSE),
   cores = 1,
@@ -168,7 +170,14 @@ Expectation Maximisation algorithm (SAEM).

 # S3 method for saem.mmkin
 print(x, digits = max(3, getOption("digits") - 3), ...)
 
-saemix_model(object, cores = 1, verbose = FALSE, ...)
+saemix_model(
+  object,
+  solution_type = "auto",
+  transformations = c("mkin", "saemix"),
+  cores = 1,
+  verbose = FALSE,
+  ...
+)
 
 saemix_data(object, verbose = FALSE, ...)

 
@@ -180,20 +189,26 @@ Expectation Maximisation algorithm (SAEM).

       
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::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 excessive forking, apparently depending on the BLAS version
-used.

+      transformations
+      
Per default, all parameter transformations are done
+in mkin. If this argument is set to 'saemix', parameter transformations
+are done in 'saemix' for the supported cases. Currently this is only
+supported in cases where the initial concentration of the parent is not fixed,
+SFO or DFOP is used for the parent and there is either no metabolite or one.

+    
+    
+      solution_type
+      
Possibility to specify the solution type in case the
+automatic choice is not desired

+    
+    
+      control
+      
Passed to saemix::saemix

     
     
       verbose
@@ -243,31 +258,31 @@ using mmkin.

 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,
+f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
   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] "Mon Nov 30 15:53:02 2020"
+#> [1] "Fri Dec 11 15:37:47 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:04 2020"

+#> [1] "Fri Dec 11 15:37:49 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] "Mon Nov 30 15:53:05 2020"
+#> [1] "Fri Dec 11 15:37:51 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:07 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#> [1] "Fri Dec 11 15:37:52 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov 30 15:53:07 2020"
+#> [1] "Fri Dec 11 15:37:52 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:09 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#> [1] "Fri Dec 11 15:37:55 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov 30 15:53:10 2020"
+#> [1] "Fri Dec 11 15:37:55 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:13 2020"

+#> [1] "Fri Dec 11 15:37:58 2020"

 # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 # functions from saemix
 library(saemix)
@@ -313,10 +328,10 @@ using mmkin.

 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] "Mon Nov 30 15:53:15 2020"
+#> [1] "Fri Dec 11 15:38:01 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:20 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+#> [1] "Fri Dec 11 15:38:06 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
 
#> Likelihoods computed by importance sampling 
#>        AIC      BIC
 #> 1 467.7096 464.9757
 #> 2 469.5208 466.3963

@@ -331,20 +346,21 @@ using mmkin.

 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
+# saem fits of SFO-SFO and DFOP-SFO to these data 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] "Mon Nov 30 15:53:23 2020"
+#> [1] "Fri Dec 11 15:38:09 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:28 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
+#> [1] "Fri Dec 11 15:38:14 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov 30 15:53:29 2020"
+#> [1] "Fri Dec 11 15:38:15 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov 30 15:53:38 2020"
# We can use print, plot and summary methods to check the results
+#> [1] "Fri Dec 11 15:38:24 2020"
# We can use print, plot and summary methods to check the results
 print(f_saem_dfop_sfo)
 
#> Kinetic nonlinear mixed-effects model fit by SAEM
 #> Structural model:
@@ -387,8 +403,8 @@ using mmkin.

 
#> 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:     Mon Nov 30 15:53:38 2020 
-#> Date of summary: Mon Nov 30 15:53:39 2020 
+#> Date of fit:     Fri Dec 11 15:38:25 2020 
+#> Date of summary: Fri Dec 11 15:38:25 2020 
 #> 
 #> Equations:
 #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -403,15 +419,15 @@ using mmkin.

 #> 
 #> Model predictions using solution type analytical 
 #> 
-#> Fitted in 9.963 s using 300, 100 iterations
+#> Fitted in 10.096 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 
+#>         93.8102         -9.7647         -0.9711         -1.8799         -4.2708 
 #>        g_qlogis 
-#>       0.1356441 
+#>          0.1356 
 #> 
 #> Fixed degradation parameter values:
 #> None
@@ -422,7 +438,7 @@ using mmkin.

 #>     AIC   BIC logLik
 #>   841.6 836.5 -407.8
 #> 
-#> Optimised, transformed parameters with symmetric confidence intervals:
+#> Optimised parameters:
 #>                     est.  lower   upper
 #> parent_0        93.76647 91.153 96.3798
 #> log_k_A1        -6.13235 -8.458 -3.8068
@@ -452,7 +468,7 @@ using mmkin.

 #>      est. lower upper
 #> a.1 1.883 1.666   2.1
 #> 
-#> Backtransformed parameters with asymmetric confidence intervals:
+#> Backtransformed parameters:
 #>                     est.     lower    upper
 #> parent_0       93.766473 9.115e+01 96.37983
 #> k_A1            0.002171 2.122e-04  0.02222
@@ -643,32 +659,19 @@ using mmkin.

 #>  Dataset 10     A1   91     13.2  11.93773 -1.26227 1.883    -0.670224
 #>  Dataset 10     A1  120     14.3  12.77666 -1.52334 1.883    -0.808847
 #>  Dataset 10     A1  120     12.1  12.77666  0.67666 1.883     0.359282

-# Using a single core, the following takes about 6 minutes as we do not have an
-# analytical solution. Using 10 cores it is slower instead of faster
-f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
-
#> Running main SAEM algorithm
-#> [1] "Mon Nov 30 15:53:39 2020"
-#> DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
-#>       taken on this call before reaching TOUT     
-#> In above message, I1 = 5000
-#>  
-#> In above message, R1 = 0.00156238
-#>  
-#> DLSODA-  At T (=R1) and step size H (=R2), the    
-#>       corrector convergence failed repeatedly     
-#>       or with ABS(H) = HMIN   
-#> In above message, R1 = 0, R2 = 1.1373e-10
-#>  
-#> DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
-#>       taken on this call before reaching TOUT     
-#> In above message, I1 = 5000
-#>  
-#> In above message, R1 = 2.24752e-06
-#>  
-#> ....
-#>     Minimisation finished
-#> [1] "Mon Nov 30 16:00:45 2020"
plot(f_saem_fomc)
-
# }
+# The following takes about 6 minutes
+#f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+#  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+
+#saemix::compare.saemix(list(
+#  f_saem_dfop_sfo$so,
+#  f_saem_dfop_sfo_deSolve$so))
+
+# If the model supports it, we can also use eigenvalue based solutions, which
+# take a similar amount of time
+#f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
+#  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+# }