Create a mixed effects model from an mmkin row object

Date: Fri, 11 Dec 2020 15:39:43 +0100 Subject: Update static docs --- _pkgdown.yml | 1 + docs/dev/news/index.html | 3 +- docs/dev/pkgdown.yml | 2 +- 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 +++++++------- docs/dev/sitemap.xml | 3 + 19 files changed, 406 insertions(+), 137 deletions(-) create mode 100644 docs/dev/reference/mixed-1.png create mode 100644 docs/dev/reference/mixed.html diff --git a/_pkgdown.yml b/_pkgdown.yml index cb3fa078..2ae14397 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -49,6 +49,7 @@ reference: - nlme_function - get_deg_func - saemix_model + - mixed - title: Datasets and known results contents: - focus_soil_moisture diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html index 70cbb00b..7059189d 100644 --- a/docs/dev/news/index.html +++ b/docs/dev/news/index.html @@ -146,7 +146,8 @@ mkin 0.9.50.4 (unreleased) Unreleased

‘mkinmod’ models gain arguments ‘name’ and ‘dll_dir’ which, in conjunction with the upcoming version of the ‘inline’ package, makes it possible to still use the DLL used for fast ODE solutions with ‘deSolve’ after saving and restoring the ‘mkinmod’ object.
‘mixed.mmkin’ New container for mmkin objects for plotting with the ‘plot.mixed.mmkin’ method
‘mkinmod’ models gain arguments ‘name’ and ‘dll_dir’ which, in conjunction with a current version of the ‘inline’ package, make it possible to still use the DLL used for fast ODE solutions with ‘deSolve’ after saving and restoring the ‘mkinmod’ object.
‘f_norm_temp_focus’ generic function to normalise time intervals using the FOCUS method, with methods for numeric vectors and ‘mkindsg’ objects
‘mkindsg’ R6 class for groups of ‘mkinds’ datasets with metadata
‘D24_2014’ dataset

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

Create a mixed effects model from an mmkin row object — mixed • mkin

+ + + + +

+ +

+ + +

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)
+
# }
+

+ +

+ + + +

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.

mkinsub

(

"SFO"

"m1"

"Test compound"

)

mkinsub

(

"SFO"

"Metabolite M1"

)

"SFO_SFO"

DLL_dir

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", ]) -