From 6464d3999338d34c081f360694dbc0bc0abf68cb Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 9 Nov 2020 14:23:16 +0100 Subject: Add plot method for saem.mmkin objects --- NAMESPACE | 2 + R/plot.saem.mmkin.R | 218 ++++++++++++++++++++++ _pkgdown.yml | 1 + docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Rplot001.png | Bin 1011 -> 27839 bytes docs/dev/reference/Rplot002.png | Bin 57363 -> 56909 bytes docs/dev/reference/index.html | 6 + docs/dev/reference/plot.nlme.mmkin.html | 10 +- docs/dev/reference/plot.saem.mmkin-1.png | Bin 0 -> 86076 bytes docs/dev/reference/plot.saem.mmkin-2.png | Bin 0 -> 164014 bytes docs/dev/reference/plot.saem.mmkin.html | 308 +++++++++++++++++++++++++++++++ docs/dev/sitemap.xml | 3 + man/plot.nlme.mmkin.Rd | 8 +- man/plot.saem.mmkin.Rd | 94 ++++++++++ 14 files changed, 637 insertions(+), 15 deletions(-) create mode 100644 R/plot.saem.mmkin.R create mode 100644 docs/dev/reference/plot.saem.mmkin-1.png create mode 100644 docs/dev/reference/plot.saem.mmkin-2.png create mode 100644 docs/dev/reference/plot.saem.mmkin.html create mode 100644 man/plot.saem.mmkin.Rd diff --git a/NAMESPACE b/NAMESPACE index ed1a46de..29462082 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -20,6 +20,7 @@ S3method(plot,mkinfit) S3method(plot,mmkin) S3method(plot,nafta) S3method(plot,nlme.mmkin) +S3method(plot,saem.mmkin) S3method(print,mkinds) S3method(print,mkinmod) S3method(print,mmkin) @@ -98,6 +99,7 @@ importFrom(parallel,mclapply) importFrom(parallel,parLapply) importFrom(pkgbuild,has_compiler) importFrom(purrr,map_dfr) +importFrom(saemix,psi) importFrom(stats,AIC) importFrom(stats,BIC) importFrom(stats,aggregate) diff --git a/R/plot.saem.mmkin.R b/R/plot.saem.mmkin.R new file mode 100644 index 00000000..ce43fdb6 --- /dev/null +++ b/R/plot.saem.mmkin.R @@ -0,0 +1,218 @@ +if(getRversion() >= '2.15.1') utils::globalVariables("ds") + +#' Plot an saem fitted nonlinear mixed model obtained via an mmkin row object +#' +#' @param x An object of class \code{\link{saem.mmkin}} +#' @param i A numeric index to select datasets for which to plot the saem fit, +#' in case plots get too large +#' @inheritParams plot.mkinfit +#' @param standardized Should the residuals be standardized? Only takes effect if +#' `resplot = "time"`. +#' @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 +#' @param ncol.legend Number of columns to use in the legend +#' @param nrow.legend Number of rows to use in the legend +#' @param resplot Should the residuals plotted against time or against +#' predicted values? +#' @param col_ds Colors used for plotting the observed data and the +#' corresponding model prediction lines for the different datasets. +#' @param pch_ds Symbols to be used for plotting the data. +#' @param lty_ds Line types to be used for the model predictions. +#' @importFrom saemix psi +#' @return 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) +#' +#' f_saem <- saem(f) +#' plot(f_saem) +#' } +#' @export +plot.saem.mmkin <- function(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, ...) +{ + + oldpar <- par(no.readonly = TRUE) + + fit_1 <- x$mmkin[[1]] + ds_names <- colnames(x$mmkin) + + degparms_optim <- psi(x$so) + rownames(degparms_optim) <- ds_names + degparms_optim_names <- setdiff(names(fit_1$par), names(fit_1$errparms)) + colnames(degparms_optim) <- degparms_optim_names + + degparms_fixed <- fit_1$fixed$value + names(degparms_fixed) <- rownames(fit_1$fixed) + degparms_all <- cbind(as.matrix(degparms_optim), + matrix(rep(degparms_fixed, nrow(degparms_optim)), + ncol = length(degparms_fixed), + nrow = nrow(degparms_optim), byrow = TRUE)) + degparms_all_names <- c(degparms_optim_names, names(degparms_fixed)) + colnames(degparms_all) <- degparms_all_names + + odeini_names <- grep("_0$", degparms_all_names, value = TRUE) + odeparms_names <- setdiff(degparms_all_names, odeini_names) + + residual_type = ifelse(standardized, "iwres", "ires") + + observed <- cbind(x$data, + residual = x$so@results@predictions[[residual_type]]) + + n_plot_rows = length(obs_vars) + n_plots = n_plot_rows * 2 + + # Set relative plot heights, so the first plot row is the norm + rel.heights <- if (n_plot_rows > 1) { + c(rel.height.legend, c(rep(1, n_plot_rows - 1), rel.height.bottom)) + } else { + c(rel.height.legend, 1) + } + + layout_matrix = matrix(c(1, 1, 2:(n_plots + 1)), + n_plot_rows + 1, 2, byrow = TRUE) + layout(layout_matrix, heights = rel.heights) + + par(mar = c(0.1, 2.1, 0.6, 2.1)) + + 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)) + + + solution_type = fit_1$solution_type + + outtimes <- sort(unique(c(x$data$time, + seq(xlim[1], xlim[2], length.out = 50)))) + + pred_ds <- purrr::map_dfr(i, function(ds_i) { + odeparms_trans <- degparms_all[ds_i, odeparms_names] + names(odeparms_trans) <- odeparms_names # needed if only one odeparm + odeparms <- backtransform_odeparms(odeparms_trans, + x$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) + + odeini <- degparms_all[ds_i, odeini_names] + names(odeini) <- gsub("_0", "", odeini_names) + + out <- mkinpredict(x$mkinmod, odeparms, odeini, + outtimes, solution_type = solution_type, + atol = fit_1$atol, rtol = fit_1$rtol) + return(cbind(as.data.frame(out), ds = ds_names[ds_i])) + }) + + degparms_pop <- x$so@results@fixed.effects + names(degparms_pop) <- degparms_optim_names + degparms_all_pop <- c(degparms_pop, degparms_fixed) + + odeparms_pop_trans <- degparms_all_pop[odeparms_names] + odeparms_pop <- backtransform_odeparms(odeparms_pop_trans, + x$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) + + odeini_pop <- degparms_all_pop[odeini_names] + names(odeini_pop) <- gsub("_0", "", odeini_names) + + pred_pop <- as.data.frame( + mkinpredict(x$mkinmod, odeparms_pop, odeini_pop, + outtimes, solution_type = solution_type, + atol = fit_1$atol, rtol = fit_1$rtol)) + + resplot <- match.arg(resplot) + + # Loop plot rows + for (plot_row in 1:n_plot_rows) { + + obs_var <- obs_vars[plot_row] + observed_row <- subset(observed, name == obs_var) + + # Set ylim to sensible default, or use ymax + if (identical(ymax, "auto")) { + ylim_row = c(0, + max(c(observed_row$value, pred_ds[[obs_var]]), na.rm = TRUE)) + } else { + ylim_row = c(0, ymax[plot_row]) + } + + # Margins for bottom row of plots when we have more than one row + # This is the only row that needs to show the x axis legend + if (plot_row == n_plot_rows) { + par(mar = c(5.1, 4.1, 2.1, 2.1)) + } else { + par(mar = c(3.0, 4.1, 2.1, 2.1)) + } + + plot(pred_pop$time, pred_pop[[obs_var]], + type = "l", lwd = 2, + xlim = xlim, ylim = ylim_row, + xlab = xlab, ylab = obs_var, frame = frame) + + 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]) + lines(subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)), + col = col_ds[ds_i], lty = lty_ds[ds_i]) + } + + if (identical(maxabs, "auto")) { + maxabs = max(abs(observed_row$residual), na.rm = TRUE) + } + + if (identical(resplot, "time")) { + plot(0, type = "n", xlim = xlim, xlab = "Time", + ylim = c(-1.2 * maxabs, 1.2 * maxabs), + ylab = if (standardized) "Standardized residual" else "Residual") + + abline(h = 0, lty = 2) + + for (ds_i in seq_along(i)) { + points(subset(observed_row, ds == ds_names[ds_i], c("time", "residual")), + col = col_ds[ds_i], pch = pch_ds[ds_i]) + } + } + + if (identical(resplot, "predicted")) { + plot(0, type = "n", + xlim = c(0, max(pred_ds[[obs_var]])), + xlab = "Predicted", + ylim = c(-1.2 * maxabs, 1.2 * maxabs), + ylab = if (standardized) "Standardized residual" else "Residual") + + abline(h = 0, lty = 2) + + for (ds_i in seq_along(i)) { + observed_row_ds <- merge( + subset(observed_row, ds == ds_names[ds_i], c("time", "residual")), + subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var))) + points(observed_row_ds[c(3, 2)], + col = col_ds[ds_i], pch = pch_ds[ds_i]) + } + } + } +} diff --git a/_pkgdown.yml b/_pkgdown.yml index b036ed59..61c69a9d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -46,6 +46,7 @@ reference: - summary.nlme.mmkin - saem.mmkin - summary.saem.mmkin + - plot.saem.mmkin - nlme_function - get_deg_func - saemix_model diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index a8bb6144..e14af5e0 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-11-09T08:23Z +last_built: 2020-11-09T13:23Z 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..cfc5bc2b 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 8ada7133..cd2014eb 100644 Binary files a/docs/dev/reference/Rplot002.png and b/docs/dev/reference/Rplot002.png differ diff --git a/docs/dev/reference/index.html b/docs/dev/reference/index.html index 825c4d27..f5621402 100644 --- a/docs/dev/reference/index.html +++ b/docs/dev/reference/index.html @@ -348,6 +348,12 @@ of an mmkin object

Summary method for class "saem.mmkin"

+ +

plot(<saem.mmkin>)

+ +

Plot an saem fitted nonlinear mixed model obtained via an mmkin row object

+ +

nlme_function() mean_degparms() nlme_data()

diff --git a/docs/dev/reference/plot.nlme.mmkin.html b/docs/dev/reference/plot.nlme.mmkin.html index 267bef05..8f5f85f7 100644 --- a/docs/dev/reference/plot.nlme.mmkin.html +++ b/docs/dev/reference/plot.nlme.mmkin.html @@ -150,7 +150,7 @@ 
# S3 method for nlme.mmkin
 plot(
   x,
-  i = 1:ncol(x$mmkin_orig),
+  i = 1:ncol(x$mmkin),
   obs_vars = names(x$mkinmod$map),
   standardized = TRUE,
   xlab = "Time",
@@ -165,8 +165,7 @@
   pch_ds = 1:length(i),
   col_ds = pch_ds + 1,
   lty_ds = col_ds,
-  frame = TRUE,
-  ...
+  frame = TRUE
 )

Arguments

@@ -247,11 +246,6 @@ corresponding model prediction lines for the different datasets.

frame

Should a frame be drawn around the plots?

- - ... -

Further arguments passed to plot.mkinfit and -mkinresplot.

-

Value

diff --git a/docs/dev/reference/plot.saem.mmkin-1.png b/docs/dev/reference/plot.saem.mmkin-1.png new file mode 100644 index 00000000..5cb33214 Binary files /dev/null and b/docs/dev/reference/plot.saem.mmkin-1.png differ diff --git a/docs/dev/reference/plot.saem.mmkin-2.png b/docs/dev/reference/plot.saem.mmkin-2.png new file mode 100644 index 00000000..67058e6c Binary files /dev/null and b/docs/dev/reference/plot.saem.mmkin-2.png differ diff --git a/docs/dev/reference/plot.saem.mmkin.html b/docs/dev/reference/plot.saem.mmkin.html new file mode 100644 index 00000000..d0790a98 --- /dev/null +++ b/docs/dev/reference/plot.saem.mmkin.html @@ -0,0 +1,308 @@ + + + + + + + + +Plot an saem fitted nonlinear mixed model obtained via an mmkin row object — plot.saem.mmkin • mkin + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Plot an saem fitted nonlinear mixed model obtained via an mmkin row object

+ Source: R/plot.saem.mmkin.R + +
+ +
+

Plot an saem fitted nonlinear mixed model obtained via an mmkin row object

+
+ + 
# S3 method for saem.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

i

A numeric index to select datasets for which to plot the saem fit, +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)
+

+f_saem <- saem(f)
+
#> Running main SAEM algorithm
+#> [1] "Mon Nov  9 14:21:22 2020"
+#> ....
+#>     Minimisation finished
+#> [1] "Mon Nov  9 14:21:31 2020"
plot(f_saem)
+
# }
+
+
+ +
+ + +
+
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 1.6.1.

+
+ +
+
+ + + + + + + + diff --git a/docs/dev/sitemap.xml b/docs/dev/sitemap.xml index d9ef2228..d6f55550 100644 --- a/docs/dev/sitemap.xml +++ b/docs/dev/sitemap.xml @@ -159,6 +159,9 @@ https://pkgdown.jrwb.de/mkin/reference/plot.nlme.mmkin.html + + https://pkgdown.jrwb.de/mkin/reference/plot.saem.mmkin.html + https://pkgdown.jrwb.de/mkin/reference/print.mkinds.html diff --git a/man/plot.nlme.mmkin.Rd b/man/plot.nlme.mmkin.Rd index 5f0f0ef1..f426f77b 100644 --- a/man/plot.nlme.mmkin.Rd +++ b/man/plot.nlme.mmkin.Rd @@ -6,7 +6,7 @@ \usage{ \method{plot}{nlme.mmkin}( x, - i = 1:ncol(x$mmkin_orig), + i = 1:ncol(x$mmkin), obs_vars = names(x$mkinmod$map), standardized = TRUE, xlab = "Time", @@ -21,8 +21,7 @@ pch_ds = 1:length(i), col_ds = pch_ds + 1, lty_ds = col_ds, - frame = TRUE, - ... + frame = TRUE ) } \arguments{ @@ -66,9 +65,6 @@ corresponding model prediction lines for the different datasets.} \item{lty_ds}{Line types to be used for the model predictions.} \item{frame}{Should a frame be drawn around the plots?} - -\item{\dots}{Further arguments passed to \code{\link{plot.mkinfit}} and -\code{\link{mkinresplot}}.} } \value{ The function is called for its side effect. diff --git a/man/plot.saem.mmkin.Rd b/man/plot.saem.mmkin.Rd new file mode 100644 index 00000000..1f674bd7 --- /dev/null +++ b/man/plot.saem.mmkin.Rd @@ -0,0 +1,94 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.saem.mmkin.R +\name{plot.saem.mmkin} +\alias{plot.saem.mmkin} +\title{Plot an saem fitted nonlinear mixed model obtained via an mmkin row object} +\usage{ +\method{plot}{saem.mmkin}( + 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{ +\item{x}{An object of class \code{\link{saem.mmkin}}} + +\item{i}{A numeric index to select datasets for which to plot the saem fit, +in case plots get too large} + +\item{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.} + +\item{standardized}{Should the residuals be standardized? Only takes effect if +\code{resplot = "time"}.} + +\item{xlab}{Label for the x axis.} + +\item{xlim}{Plot range in x direction.} + +\item{resplot}{Should the residuals plotted against time or against +predicted values?} + +\item{ymax}{Vector of maximum y axis values} + +\item{maxabs}{Maximum absolute value of the residuals. This is used for the +scaling of the y axis and defaults to "auto".} + +\item{ncol.legend}{Number of columns to use in the legend} + +\item{nrow.legend}{Number of rows to use in the legend} + +\item{rel.height.legend}{The relative height of the legend shown on top} + +\item{rel.height.bottom}{The relative height of the bottom plot row} + +\item{pch_ds}{Symbols to be used for plotting the data.} + +\item{col_ds}{Colors used for plotting the observed data and the +corresponding model prediction lines for the different datasets.} + +\item{lty_ds}{Line types to be used for the model predictions.} + +\item{frame}{Should a frame be drawn around the plots?} + +\item{...}{Further arguments passed to \code{\link{plot}}.} +} +\value{ +The function is called for its side effect. +} +\description{ +Plot an saem fitted nonlinear mixed model obtained via an mmkin row object +} +\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) + +f_saem <- saem(f) +plot(f_saem) +} +} +\author{ +Johannes Ranke +} -- cgit v1.2.1