Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object

# S3 method for mmkin -nlme( +nlme( model, data = sys.frame(sys.parent()), fixed, @@ -253,8 +253,8 @@ parameters taken from the mmkin object are used

Value

Upon success, a fitted nlme.mmkin object, which is an nlme object -with additional elements

Upon success, a fitted 'nlme.mmkin' object, which is an nlme object +with additional elements. It also inherits from 'mixed.mmkin'.

Note

As the object inherits from nlme::nlme, there is a wealth of @@ -262,7 +262,8 @@ methods that will automatically work on 'nlme.mmkin' objects, such as nlme::intervals(), nlme::anova.lme() and nlme::coef.lme().

Examples

ds <- lapply(experimental_data_for_UBA_2019[6:10],
@@ -288,21 +289,21 @@ methods that will automatically work on 'nlme.mmkin' objects, such as
 #> Log-likelihood: -238.5635
 #> 
 #> Fixed effects:
-#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
-#>    parent_0      log_k1      log_k2       g_ilr 
-#> 94.17015133 -1.80015306 -4.14738870  0.02290935 
+#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) 
+#>    parent_0      log_k1      log_k2    g_qlogis 
+#> 94.17015185 -1.80015278 -4.14738834  0.03239833 
 #> 
 #> Random effects:
-#>  Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1)
+#>  Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1)
 #>  Level: ds
 #>  Structure: Diagonal
-#>         parent_0    log_k1  log_k2     g_ilr Residual
-#> StdDev: 2.488249 0.8447273 1.32965 0.3289311 2.321364
+#>         parent_0    log_k1  log_k2  g_qlogis Residual
+#> StdDev: 2.488249 0.8447275 1.32965 0.4651789 2.321364
 #> 
plot(f_nlme_dfop)
 
endpoints(f_nlme_dfop)
 
#> $distimes
 #>            DT50     DT90 DT50back  DT50_k1  DT50_k2
-#> parent 10.79857 100.7937 30.34193 4.193938 43.85443
+#> parent 10.79857 100.7937 30.34192 4.193937 43.85442
 #> 

 # \dontrun{
   f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1))
@@ -353,35 +354,35 @@ methods that will automatically work on 'nlme.mmkin' objects, such as
     control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE)
 
#> 
 #> **Iteration 1
-#> LME step: Loglik: -404.9582, nlminb iterations: 1
+#> LME step: Loglik: -404.9583, nlminb iterations: 1
 #> reStruct  parameters:
 #>        ds1        ds2        ds3        ds4        ds5        ds6 
-#> -0.4114355  0.9798697  1.6990037  0.7293315  0.3354323  1.7113046 
+#> -0.4114356  0.9798646  1.3524300  0.7293315  0.3354323  1.3647313 
 #>  Beginning PNLS step: ..  completed fit_nlme() step.
-#> PNLS step: RSS =  630.3644 
-#>  fixed effects: 93.82269  -5.455991  -0.6788957  -1.862196  -4.199671  0.05532828  
+#> PNLS step: RSS =  630.3633 
+#>  fixed effects: 93.82269  -5.455993  -0.9601037  -1.862196  -4.199671  0.07824609  
 #>  iterations: 120 
 #> Convergence crit. (must all become <= tolerance = 0.0005):
 #>     fixed  reStruct 
-#> 0.7885368 0.5822683 
+#> 0.7897284 0.5822782 
 #> 
 #> **Iteration 2
 #> LME step: Loglik: -407.7755, nlminb iterations: 11
 #> reStruct  parameters:
-#>          ds1          ds2          ds3          ds4          ds5          ds6 
-#> -0.371224133  0.003056179  1.789939402  0.724671158  0.301602977  1.754200729 
+#>         ds1         ds2         ds3         ds4         ds5         ds6 
+#> -0.37122411  0.00305562  1.44336560  0.72467122  0.30160310  1.40762692 
 #>  Beginning PNLS step: ..  completed fit_nlme() step.
-#> PNLS step: RSS =  630.3633 
-#>  fixed effects: 93.82269  -5.455992  -0.6788958  -1.862196  -4.199671  0.05532831  
+#> PNLS step: RSS =  630.3637 
+#>  fixed effects: 93.82269  -5.455992  -0.9601036  -1.862196  -4.199671  0.0782462  
 #>  iterations: 120 
 #> Convergence crit. (must all become <= tolerance = 0.0005):
 #>        fixed     reStruct 
-#> 4.789774e-07 2.200661e-05 

+#> 1.375673e-06 9.758294e-06 

   plot(f_nlme_dfop_sfo)
 

   anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo)
 
#>                 Model df       AIC       BIC    logLik   Test  L.Ratio p-value
-#> f_nlme_dfop_sfo     1 13  843.8547  884.6201 -408.9273                        
+#> f_nlme_dfop_sfo     1 13  843.8547  884.6201 -408.9274                        
 #> f_nlme_sfo_sfo      2  9 1085.1821 1113.4043 -533.5910 1 vs 2 249.3274  <.0001

   endpoints(f_nlme_sfo_sfo)
 
#> $ff
@@ -400,7 +401,7 @@ methods that will automatically work on 'nlme.mmkin' objects, such as
 #> $distimes
 #>             DT50     DT90 DT50back  DT50_k1  DT50_k2
 #> parent  11.07091 104.6320 31.49738 4.462384 46.20825
-#> A1     162.30536 539.1667       NA       NA       NA
+#> A1     162.30523 539.1663       NA       NA       NA
 #> 

   if (length(findFunction("varConstProp")) > 0) { # tc error model for nlme available
     # Attempts to fit metabolite kinetics with the tc error model are possible,
@@ -425,23 +426,23 @@ methods that will automatically work on 'nlme.mmkin' objects, such as
 #> Log-likelihood: -238.4298
 #> 
 #> Fixed effects:
-#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
-#>    parent_0      log_k1      log_k2       g_ilr 
-#> 94.04774463 -1.82339924 -4.16715509  0.04020161 
+#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) 
+#>    parent_0      log_k1      log_k2    g_qlogis 
+#> 94.04774566 -1.82339808 -4.16715311  0.05685186 
 #> 
 #> Random effects:
-#>  Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1)
+#>  Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1)
 #>  Level: ds
 #>  Structure: Diagonal
-#>         parent_0    log_k1   log_k2     g_ilr Residual
-#> StdDev: 2.473883 0.8499901 1.337187 0.3294411        1
+#>         parent_0    log_k1   log_k2  g_qlogis Residual
+#> StdDev: 2.473881 0.8499884 1.337185 0.4659005        1
 #> 
 #> Variance function:
 #>  Structure: Constant plus proportion of variance covariate
 #>  Formula: ~fitted(.) 
 #>  Parameter estimates:
 #>      const       prop 
-#> 2.23222625 0.01262414 

+#> 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")
diff --git a/docs/dev/reference/plot.mixed.mmkin-1.png b/docs/dev/reference/plot.mixed.mmkin-1.png
new file mode 100644
index 00000000..5cb33214
Binary files /dev/null and b/docs/dev/reference/plot.mixed.mmkin-1.png differ
diff --git a/docs/dev/reference/plot.mixed.mmkin-2.png b/docs/dev/reference/plot.mixed.mmkin-2.png
new file mode 100644
index 00000000..c0d67204
Binary files /dev/null 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
new file mode 100644
index 00000000..67058e6c
Binary files /dev/null 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
new file mode 100644
index 00000000..4108aea3
--- /dev/null
+++ b/docs/dev/reference/plot.mixed.mmkin.html
@@ -0,0 +1,313 @@
+
+
+
+  
+  
+
+
+
+Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object — plot.mixed.mmkin • mkin
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+  
+
+  
+    
+      
+      
+  
+    
+      
+      
+        mkin
+        0.9.50.4
+      
+    
+
+    
+      
+        
+  Functions and data
+
+
+  
+    Articles
+     
+    
+  
+  
+    
+      Introduction to mkin
+    
+    
+      Example evaluation of FOCUS Example Dataset D
+    
+    
+      Example evaluation of FOCUS Laboratory Data L1 to L3
+    
+    
+      Example evaluation of FOCUS Example Dataset Z
+    
+    
+      Performance benefit by using compiled model definitions in mkin
+    
+    
+      Calculation of time weighted average concentrations with mkin
+    
+    
+      Example evaluation of NAFTA SOP Attachment examples
+    
+    
+      Some benchmark timings
+    
+  
+
+
+  News
+
+      
+      
+        
+  
+    
+     
+  
+
+      
+      
+    

+  

+

+
+      
+
+      
+
+
+  
+    
+    Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
+    Source: R/plot.mixed.mmkin.R
+    plot.mixed.mmkin.Rd
+    
+
+    
+    Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
+    
+
+    # S3 method for mixed.mmkin
+plot(
+  x,
+  i = 1:ncol(x$mmkin),
+  obs_vars = names(x$mkinmod$map),
+  standardized = TRUE,
+  xlab = "Time",
+  xlim = range(x$data$time),
+  resplot = c("predicted", "time"),
+  ymax = "auto",
+  maxabs = "auto",
+  ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
+  nrow.legend = ceiling((length(i) + 1)/ncol.legend),
+  rel.height.legend = 0.03 + 0.08 * nrow.legend,
+  rel.height.bottom = 1.1,
+  pch_ds = 1:length(i),
+  col_ds = pch_ds + 1,
+  lty_ds = col_ds,
+  frame = TRUE,
+  ...
+)
+
+    Arguments
+    
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    
+      
+      
+    
+    x An object of class saem.mmkin or nlme.mmkin
i A numeric index to select datasets for which to plot the individual predictions,
+in case plots get too large
obs_vars A character vector of names of the observed variables for
+which the data and the model should be plotted. Defauls to all observed
+variables in the model.
standardized Should the residuals be standardized? Only takes effect if
+resplot = "time".
xlab Label for the x axis.
xlim Plot range in x direction.
resplot Should the residuals plotted against time or against
+predicted values?
ymax Vector of maximum y axis values
maxabs Maximum absolute value of the residuals. This is used for the
+scaling of the y axis and defaults to "auto".
ncol.legend Number of columns to use in the legend
nrow.legend Number of rows to use in the legend
rel.height.legend The relative height of the legend shown on top
rel.height.bottom The relative height of the bottom plot row
pch_ds Symbols to be used for plotting the data.
col_ds Colors used for plotting the observed data and the
+corresponding model prediction lines for the different datasets.
lty_ds Line types to be used for the model predictions.
frame Should a frame be drawn around the plots?
... Further arguments passed to plot.
+
+    Value
+
+    The function is called for its side effect.
+    Author
+
+    Johannes Ranke
+
+    Examples
+    ds <- lapply(experimental_data_for_UBA_2019[6:10],
+ function(x) x$data[c("name", "time", "value")])
+names(ds) <- paste0("ds ", 6:10)
+dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+  A1 = mkinsub("SFO"), quiet = TRUE)
+# \dontrun{
+f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
+plot(f[, 3:4], standardized = TRUE)
+

+# 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))
+plot(f_nlme)
+

+f_saem <- saem(f)
+
#> Running main SAEM algorithm
+#> [1] "Mon Nov  9 17:18:17 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Mon Nov  9 17:18:26 2020"
plot(f_saem)
+
# }
+
+  
+  
+    
+      Contents
+    
+  
+
+
+
+      
+      
+  Developed by Johannes Ranke.
+
+
+
+  Site built with pkgdown 1.6.1.
+
+
+      
+   
+
+  
+
+
+  
+
+
+
diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html
index 26f4c3e3..f1b4c421 100644
--- a/docs/dev/reference/saem.html
+++ b/docs/dev/reference/saem.html
@@ -206,7 +206,8 @@ by the saemix function?
     Value
 
     An S3 object of class 'saem.mmkin', containing the fitted
-saemix::SaemixObject as a list component named 'so'.
+saemix::SaemixObject as a list component named 'so'. The
+object also inherits from 'mixed.mmkin'.
 An saemix::SaemixModel object.
 An saemix::SaemixData object.
     Details
@@ -218,7 +219,7 @@ psi0 of saemix::
 using mmkin.
     See also
 
-    summary.saem.mmkin
+    summary.saem.mmkin plot.mixed.mmkin
 
     Examples
     # \dontrun{
@@ -229,27 +230,27 @@ using mmkin.
   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  9 09:03:11 2020"
+#> [1] "Mon Nov  9 17:18:28 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:13 2020"

+#> [1] "Mon Nov  9 17:18:30 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  9 09:03:14 2020"
+#> [1] "Mon Nov  9 17:18:31 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:16 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#> [1] "Mon Nov  9 17:18:33 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov  9 09:03:16 2020"
+#> [1] "Mon Nov  9 17:18:33 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:18 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#> [1] "Mon Nov  9 17:18:35 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov  9 09:03:19 2020"
+#> [1] "Mon Nov  9 17:18:36 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:22 2020"

+#> [1] "Mon Nov  9 17:18:39 2020"

 # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 # functions from saemix
 library(saemix)
@@ -295,10 +296,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  9 09:03:24 2020"
+#> [1] "Mon Nov  9 17:18:41 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:29 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+#> [1] "Mon Nov  9 17:18:46 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

@@ -318,20 +319,20 @@ using mmkin.
 # 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  9 09:03:31 2020"
+#> [1] "Mon Nov  9 17:18:48 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:36 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
+#> [1] "Mon Nov  9 17:18:53 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
 
#> Running main SAEM algorithm
-#> [1] "Mon Nov  9 09:03:37 2020"
+#> [1] "Mon Nov  9 17:18:54 2020"
 #> ....
 #>     Minimisation finished
-#> [1] "Mon Nov  9 09:03:46 2020"
summary(f_saem_dfop_sfo, data = FALSE)
+#> [1] "Mon Nov  9 17:19:03 2020"
summary(f_saem_dfop_sfo, data = FALSE)
 
#> 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  9 09:03:47 2020 
-#> Date of summary: Mon Nov  9 09:03:47 2020 
+#> Date of fit:     Mon Nov  9 17:19:04 2020 
+#> Date of summary: Mon Nov  9 17:19:04 2020 
 #> 
 #> Equations:
 #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -346,7 +347,7 @@ using mmkin.
 #> 
 #> Model predictions using solution type analytical 
 #> 
-#> Fitted in 9.758 s using 300, 100 iterations
+#> Fitted in 9.941 s using 300, 100 iterations
 #> 
 #> Variance model: Constant variance 
 #> 
-- 
cgit v1.2.1

x	An object of class saem.mmkin or nlme.mmkin
i	A numeric index to select datasets for which to plot the individual predictions, +in case plots get too large
obs_vars	A character vector of names of the observed variables for +which the data and the model should be plotted. Defauls to all observed +variables in the model.
standardized	Should the residuals be standardized? Only takes effect if +`resplot = "time"`.
xlab	Label for the x axis.
xlim	Plot range in x direction.
resplot	Should the residuals plotted against time or against +predicted values?
ymax	Vector of maximum y axis values
maxabs	Maximum absolute value of the residuals. This is used for the +scaling of the y axis and defaults to "auto".
ncol.legend	Number of columns to use in the legend
nrow.legend	Number of rows to use in the legend
rel.height.legend	The relative height of the legend shown on top
rel.height.bottom	The relative height of the bottom plot row
pch_ds	Symbols to be used for plotting the data.
col_ds	Colors used for plotting the observed data and the +corresponding model prediction lines for the different datasets.
lty_ds	Line types to be used for the model predictions.
frame	Should a frame be drawn around the plots?
...	Further arguments passed to `plot`.