From d28ce9f8ad6f9573e403ebd8eb637ecd5e5b0e02 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Mon, 21 Dec 2020 06:02:10 +0100
Subject: plot.mixed: Possibility to overlay predictions
---
R/plot.mixed.mmkin.R | 37 +++++++++++++++----
docs/dev/pkgdown.yml | 2 +-
docs/dev/reference/Rplot001.png | Bin 1011 -> 19395 bytes
docs/dev/reference/plot.mixed.mmkin-3.png | Bin 163488 -> 163536 bytes
docs/dev/reference/plot.mixed.mmkin-4.png | Bin 0 -> 166687 bytes
docs/dev/reference/plot.mixed.mmkin.html | 23 +++++++++---
docs/dev/reference/saem.html | 57 +++++++++++++++---------------
man/plot.mixed.mmkin.Rd | 13 ++++++-
man/saem.Rd | 5 +++
9 files changed, 95 insertions(+), 42 deletions(-)
create mode 100644 docs/dev/reference/plot.mixed.mmkin-4.png
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index db29376e..109df283 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -8,6 +8,8 @@ utils::globalVariables("ds")
#' @inheritParams plot.mkinfit
#' @param standardized Should the residuals be standardized? Only takes effect if
#' `resplot = "time"`.
+#' @param pred_over Named list of alternative predictions as obtained
+#' from [mkinpredict] with a compatible [mkinmod].
#' @param rel.height.legend The relative height of the legend shown on top
#' @param rel.height.bottom The relative height of the bottom plot row
#' @param ymax Vector of maximum y axis values
@@ -37,8 +39,15 @@ utils::globalVariables("ds")
#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
#' plot(f_nlme)
#'
-#' f_saem <- saem(f)
+#' f_saem <- saem(f, transformations = "saemix")
#' plot(f_saem)
+#'
+#' # We can overlay the two variants if we generate predictions
+#' pred_nlme <- mkinpredict(dfop_sfo,
+#' f_nlme$bparms.optim[-1],
+#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
+#' seq(0, 180, by = 0.2))
+#' plot(f_saem, pred_over = list(nlme = pred_nlme))
#' }
#' @export
plot.mixed.mmkin <- function(x,
@@ -48,6 +57,7 @@ plot.mixed.mmkin <- function(x,
xlab = "Time",
xlim = range(x$data$time),
resplot = c("predicted", "time"),
+ pred_over = NULL,
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),
@@ -174,12 +184,19 @@ plot.mixed.mmkin <- function(x,
par(mar = c(0.1, 2.1, 0.6, 2.1))
+ # Empty plot with legend
+ if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
+ else lty_over <- NULL
+ n_pop <- 1 + length(lty_over)
+ lty_pop <- c(1, lty_over)
+
plot(0, type = "n", axes = FALSE, ann = FALSE)
legend("center", bty = "n", ncol = ncol.legend,
- legend = c("Population", ds_names[i]),
- lty = c(1, lty_ds), lwd = c(2, rep(1, length(i))),
- col = c(1, col_ds),
- pch = c(NA, pch_ds))
+ legend = c("Population", names(pred_over), ds_names[i]),
+ lty = c(lty_pop, lty_ds),
+ lwd = c(rep(2, n_pop), rep(1, length(i))),
+ col = c(rep(1, n_pop), col_ds),
+ pch = c(rep(NA, n_pop), pch_ds))
resplot <- match.arg(resplot)
@@ -206,10 +223,18 @@ plot.mixed.mmkin <- function(x,
}
plot(pred_pop$time, pred_pop[[obs_var]],
- type = "l", lwd = 2,
+ type = "l", lwd = 2, lty = lty_pop,
xlim = xlim, ylim = ylim_row,
xlab = xlab, ylab = obs_var, frame = frame)
+ if (!is.null(pred_over)) {
+ for (i_over in seq_along(pred_over)) {
+ pred_frame <- as.data.frame(pred_over[[i_over]])
+ lines(pred_frame$time, pred_frame[[obs_var]],
+ lwd = 2, lty = lty_over[i_over])
+ }
+ }
+
for (ds_i in seq_along(i)) {
points(subset(observed_row, ds == ds_names[ds_i], c("time", "value")),
col = col_ds[ds_i], pch = pch_ds[ds_i])
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index b3124589..85a358df 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -10,7 +10,7 @@ articles:
web_only/NAFTA_examples: NAFTA_examples.html
web_only/benchmarks: benchmarks.html
web_only/compiled_models: compiled_models.html
-last_built: 2020-12-19T11:24Z
+last_built: 2020-12-21T04:57Z
urls:
reference: https://pkgdown.jrwb.de/mkin/reference
article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png
index 17a35806..bca41e2c 100644
Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ
diff --git a/docs/dev/reference/plot.mixed.mmkin-3.png b/docs/dev/reference/plot.mixed.mmkin-3.png
index 3055c0c9..5e00afe6 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-4.png b/docs/dev/reference/plot.mixed.mmkin-4.png
new file mode 100644
index 00000000..6a5f3b9c
Binary files /dev/null and b/docs/dev/reference/plot.mixed.mmkin-4.png differ
diff --git a/docs/dev/reference/plot.mixed.mmkin.html b/docs/dev/reference/plot.mixed.mmkin.html
index 90ba9184..55c411e7 100644
--- a/docs/dev/reference/plot.mixed.mmkin.html
+++ b/docs/dev/reference/plot.mixed.mmkin.html
@@ -121,7 +121,7 @@
-
-
+
@@ -156,6 +156,7 @@
xlab = "Time",
xlim = range(x$data$time),
resplot = c("predicted", "time"),
+ pred_over = NULL,
ymax = "auto",
maxabs = "auto",
ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
@@ -204,6 +205,11 @@ variables in the model.
resplot |
Should the residuals plotted against time or against
predicted values? |
+
+
+ pred_over |
+ Named list of alternative predictions as obtained
+from mkinpredict with a compatible mkinmod. |
ymax |
@@ -275,13 +281,20 @@ corresponding model prediction lines for the different datasets.
f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
plot(f_nlme)
#> Running main SAEM algorithm
-#> [1] "Fri Dec 11 15:37:36 2020"
+#> [1] "Mon Dec 21 05:58:23 2020"
#> ....
#> Minimisation finished
-#> [1] "Fri Dec 11 15:37:45 2020"
# }
+#> [1] "Mon Dec 21 05:58:30 2020"
# }
+
+ degparms_start |
+ Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object. |
solution_type |
@@ -260,36 +267,33 @@ 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] "Wed Dec 16 18:27:30 2020"
+#> [1] "Mon Dec 21 05:58:33 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:31 2020"
+#> [1] "Mon Dec 21 05:58:34 2020"
#> Running main SAEM algorithm
-#> [1] "Wed Dec 16 18:27:33 2020"
+#> [1] "Mon Dec 21 05:58:36 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:35 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#> [1] "Mon Dec 21 05:58:37 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
#> Running main SAEM algorithm
-#> [1] "Wed Dec 16 18:27:35 2020"
+#> [1] "Mon Dec 21 05:58:38 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:37 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#> [1] "Mon Dec 21 05:58:40 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
#> Running main SAEM algorithm
-#> [1] "Wed Dec 16 18:27:37 2020"
+#> [1] "Mon Dec 21 05:58:40 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:40 2020"
+#> [1] "Mon Dec 21 05:58:43 2020"
#> Package saemix, version 3.1.9000
#> please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
#> Likelihoods computed by importance sampling
#> AIC BIC
-#> 1 624.2484 622.2956
-#> 2 467.7096 464.9757
-#> 3 495.4373 491.9222
#> Error in compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)): 'compare.saemix' requires at least two models.
#> Plotting convergence plots
#> Plotting individual fits
#> Simulating data using nsim = 1000 simulated datasets
@@ -326,13 +330,11 @@ 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] "Wed Dec 16 18:27:43 2020"
+#> [1] "Mon Dec 21 05:58:46 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:48 2020"
#> Likelihoods computed by importance sampling
#> AIC BIC
-#> 1 467.7096 464.9757
-#> 2 469.5208 466.3963
+#> [1] "Mon Dec 21 05:58:51 2020"
#> Error in compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)): 'compare.saemix' requires at least two models.
#> Temporary DLL for differentials generated and loaded
#> Running main SAEM algorithm
-#> [1] "Wed Dec 16 18:27:50 2020"
+#> [1] "Mon Dec 21 05:58:53 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:27:55 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
+#> [1] "Mon Dec 21 05:58:58 2020"
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
#> Running main SAEM algorithm
-#> [1] "Wed Dec 16 18:27:56 2020"
+#> [1] "Mon Dec 21 05:58:59 2020"
#> ....
#> Minimisation finished
-#> [1] "Wed Dec 16 18:28:05 2020"
# We can use print, plot and summary methods to check the results
+#> [1] "Mon Dec 21 05:59:08 2020"
#> Kinetic nonlinear mixed-effects model fit by SAEM
#> Structural model:
@@ -373,8 +375,6 @@ using
mmkin.
#> 170 observations of 2 variable(s) grouped in 5 datasets
#>
#> Likelihood computed by importance sampling
-#>
-#> LL by is "-407.78 (df=13)"
#> AIC BIC logLik
#> 841.6 836.5 -407.8
#>
@@ -400,12 +400,11 @@ using
mmkin.
#> SD.log_k2 1.90634 0.70934 3.1033
#> SD.g_qlogis 0.44771 -0.86417 1.7596
#>
-#> LL by is "-407.78 (df=13)"
#> saemix version used for fitting: 3.1.9000
+
#> 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: Wed Dec 16 18:28:06 2020
-#> Date of summary: Wed Dec 16 18:28:06 2020
+#> Date of fit: Mon Dec 21 05:59:09 2020
+#> Date of summary: Mon Dec 21 05:59:09 2020
#>
#> Equations:
#> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -420,7 +419,7 @@ using
mmkin.
#>
#> Model predictions using solution type analytical
#>
-#> Fitted in 9.889 s using 300, 100 iterations
+#> Fitted in 9.836 s using 300, 100 iterations
#>
#> Variance model: Constant variance
#>
diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd
index b90a4b3a..c7b2344f 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -12,6 +12,7 @@
xlab = "Time",
xlim = range(x$data$time),
resplot = c("predicted", "time"),
+ pred_over = NULL,
ymax = "auto",
maxabs = "auto",
ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
@@ -45,6 +46,9 @@ variables in the model.}
\item{resplot}{Should the residuals plotted against time or against
predicted values?}
+\item{pred_over}{Named list of alternative predictions as obtained
+from \link{mkinpredict} with a compatible \link{mkinmod}.}
+
\item{ymax}{Vector of maximum y axis values}
\item{maxabs}{Maximum absolute value of the residuals. This is used for the
@@ -90,8 +94,15 @@ plot(f[, 3:4], standardized = TRUE)
f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
plot(f_nlme)
-f_saem <- saem(f)
+f_saem <- saem(f, transformations = "saemix")
plot(f_saem)
+
+# We can overlay the two variants if we generate predictions
+pred_nlme <- mkinpredict(dfop_sfo,
+ f_nlme$bparms.optim[-1],
+ c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
+ seq(0, 180, by = 0.2))
+plot(f_saem, pred_over = list(nlme = pred_nlme))
}
}
\author{
diff --git a/man/saem.Rd b/man/saem.Rd
index dba948bb..775a8e7b 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -13,6 +13,7 @@ saem(object, ...)
\method{saem}{mmkin}(
object,
transformations = c("mkin", "saemix"),
+ degparms_start = numeric(),
solution_type = "auto",
control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs =
FALSE),
@@ -28,6 +29,7 @@ saemix_model(
object,
solution_type = "auto",
transformations = c("mkin", "saemix"),
+ degparms_start = numeric(),
verbose = FALSE,
...
)
@@ -46,6 +48,9 @@ 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.}
+\item{degparms_start}{Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object.}
+
\item{solution_type}{Possibility to specify the solution type in case the
automatic choice is not desired}
--
cgit v1.2.1