From dd80ab5d64bc6b56587b0542dc95e965a3a25590 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 24 Oct 2020 00:32:00 +0200 Subject: Plot method for nlme.mmkin objects Update docs --- NAMESPACE | 1 + NEWS.md | 2 + R/nlme.mmkin.R | 2 + R/plot.mmkin.R | 18 ++- R/plot.nlme.mmkin.R | 259 ++++++++++++++++++++++--------- _pkgdown.yml | 2 + docs/dev/news/index.html | 4 +- docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Extract.mmkin.html | 71 +++++---- docs/dev/reference/Rplot001.png | Bin 25568 -> 27839 bytes docs/dev/reference/Rplot002.png | Bin 26000 -> 63682 bytes docs/dev/reference/index.html | 12 ++ docs/dev/reference/plot.nlme.mmkin-1.png | Bin 35382 -> 86071 bytes docs/dev/reference/plot.nlme.mmkin-2.png | Bin 35346 -> 177237 bytes docs/dev/reference/plot.nlme.mmkin.html | 141 +++++++++++------ docs/dev/reference/print.mmkin.html | 194 +++++++++++++++++++++++ docs/dev/sitemap.xml | 3 + man/Extract.mmkin.Rd | 3 +- man/plot.nlme.mmkin.Rd | 83 ++++++---- man/print.mmkin.Rd | 16 ++ 20 files changed, 624 insertions(+), 189 deletions(-) create mode 100644 docs/dev/reference/print.mmkin.html create mode 100644 man/print.mmkin.Rd diff --git a/NAMESPACE b/NAMESPACE index 0a04d6c3..645c145e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,7 @@ S3method(plot,nafta) S3method(plot,nlme.mmkin) S3method(print,mkinds) S3method(print,mkinmod) +S3method(print,mmkin) S3method(print,nafta) S3method(print,nlme.mmkin) S3method(print,summary.mkinfit) diff --git a/NEWS.md b/NEWS.md index b7e8f38c..ea4b5318 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # mkin 0.9.50.4 (unreleased) +- 'plot' method for 'nlme.mmkin' objects + - 'print' method for 'mmkin' objects - 'saemix_model', 'saemix_data': Helper functions to fit nonlinear mixed-effects models for mmkin row objects using the saemix package diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index c8a99d59..d4720e79 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -196,6 +196,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), val <- do.call("nlme.formula", thisCall) val$mmkin_orig <- model + val$data <- thisCall[["data"]] + val$mkinmod <- model[[1]]$mkinmod class(val) <- c("nlme.mmkin", "nlme", "lme") return(val) } diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R index ac66e542..0523e4d3 100644 --- a/R/plot.mmkin.R +++ b/R/plot.mmkin.R @@ -54,12 +54,16 @@ #' #' @export plot.mmkin <- function(x, main = "auto", legends = 1, - resplot = c("time", "errmod"), - standardized = FALSE, - show_errmin = TRUE, - errmin_var = "All data", errmin_digits = 3, - cex = 0.7, rel.height.middle = 0.9, - ymax = "auto", ...) { + resplot = c("time", "errmod"), + standardized = FALSE, + show_errmin = TRUE, + errmin_var = "All data", errmin_digits = 3, + cex = 0.7, rel.height.middle = 0.9, + ymax = "auto", ...) +{ + + oldpar <- par(no.readonly = TRUE) + n.m <- nrow(x) n.d <- ncol(x) @@ -82,8 +86,6 @@ plot.mmkin <- function(x, main = "auto", legends = 1, datasets = rownames(x)) } - oldpar <- par(no.readonly = TRUE) - # Set relative plot heights, so the first and the last plot are the norm # and the middle plots (if n.fits >2) are smaller by rel.height.middle rel.heights <- if (n.fits > 2) c(1, rep(rel.height.middle, n.fits - 2), 1) diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R index 0f3ad715..084099ac 100644 --- a/R/plot.nlme.mmkin.R +++ b/R/plot.nlme.mmkin.R @@ -1,115 +1,224 @@ +if(getRversion() >= '2.15.1') utils::globalVariables("ds") + #' Plot a fitted nonlinear mixed model obtained via an mmkin row object #' #' @param x An object of class \code{\link{nlme.mmkin}} #' @param i A numeric index to select datasets for which to plot the nlme fit, #' in case plots get too large #' @param main The main title placed on the outer margin of the plot. +#' @inheritParams plot.mkinfit #' @param legends An index for the fits for which legends should be shown. -#' @param resplot Should the residuals plotted against time, using -#' \code{\link{mkinresplot}}, or as squared residuals against predicted -#' values, with the error model, using \code{\link{mkinerrplot}}. -#' @param standardized Should the residuals be standardized? This option -#' is passed to \code{\link{mkinresplot}}, it only takes effect if +#' @param standardized Should the residuals be standardized? Only takes effect if #' `resplot = "time"`. -#' @param show_errmin Should the chi2 error level be shown on top of the plots -#' to the left? -#' @param errmin_var The variable for which the FOCUS chi2 error value should -#' be shown. -#' @param errmin_digits The number of significant digits for rounding the FOCUS -#' chi2 error percentage. #' @param cex Passed to the plot functions and \code{\link{mtext}}. #' @param rel.height.middle The relative height of the middle plot, if more #' than two rows of plots are shown. -#' @param ymax Maximum y axis value for \code{\link{plot.mkinfit}}. +#' @param ymax Vector of maximum y axis values #' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and #' \code{\link{mkinresplot}}. +#' @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 stats coefficients #' @return The function is called for its side effect. #' @author Johannes Ranke #' @examples #' ds <- lapply(experimental_data_for_UBA_2019[6:10], -#' function(x) subset(x$data[c("name", "time", "value")], name == "parent")) -#' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) -#' #plot(f) # too many panels for pkgdown -#' plot(f[, 3:4]) +#' 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) +#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1) +#' plot(f[, 3:4], standardized = TRUE) #' library(nlme) -#' f_nlme <- nlme(f) -#' -#' #plot(f_nlme) # too many panels for pkgdown -#' plot(f_nlme, 3:4) +#' # 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-4)) +#' plot(f_nlme) #' @export plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), - main = "auto", legends = 1, - resplot = c("time", "errmod"), - standardized = FALSE, - show_errmin = TRUE, - errmin_var = "All data", errmin_digits = 3, + main = rownames(x$mmkin_orig), + obs_vars = names(x$mkinmod$map), + standardized = TRUE, + xlab = "Time", ylab = "Observed", + xlim = range(x$data$time), + legends = 1, + lpos = "topright", inset = c(0.05, 0.05), + resplot = c("predicted", "time"), + ymax = "auto", maxabs = "auto", cex = 0.7, rel.height.middle = 0.9, - ymax = "auto", ...) + pch_ds = 1:length(i), + col_ds = pch_ds + 1, + lty_ds = col_ds, + frame = TRUE, ...) { - degparms_optim_nlme <- coefficients(x) - degparms_optim_names <- names(degparms_optim_nlme) + oldpar <- par(no.readonly = TRUE) - odeini_optim_names <- grep("_0$", degparms_optim_names, value = TRUE) - odeparms_optim_names <- setdiff(degparms_optim_names, odeini_optim_names) + fit_1 = x$mmkin_orig[[1]] + ds_names <- colnames(x$mmkin_orig) - fit_1 <- x$mmkin_orig[[1]] + degparms_optim <- coefficients(x) + degparms_optim_names <- names(degparms_optim) + degparms_fixed <- fit_1$fixed$value + names(degparms_fixed) <- rownames(fit_1$fixed) + degparms_all <- cbind(as.matrix(degparms_optim), fit_1$bparms.fixed) + degparms_all_names <- c(degparms_optim_names, names(degparms_fixed)) + colnames(degparms_all) <- degparms_all_names - mkinfit_call <- as.list(fit_1$call)[-1] - mkinfit_call[["mkinmod"]] <- fit_1$mkinmod + odeini_names <- grep("_0$", degparms_all_names, value = TRUE) + odeparms_names <- setdiff(degparms_all_names, odeini_names) - ds <- lapply(x$mmkin_orig, function(x) { - data.frame(name = x$data$variable, - time = x$data$time, - value = x$data$observed) - }) + residual_type = ifelse(standardized, "pearson", "response") - # This takes quite some time. This could be greatly reduced - # if the plot.mkinfit code would be imported and adapted, - # allowing also to overly plots of mmkin fits and nlme fits - mmkin_nlme <- lapply(i, function(a) { + observed <- cbind(x$data, + residual = residuals(x, type = residual_type)) - degparms_optim <- as.numeric(degparms_optim_nlme[a, ]) - names(degparms_optim) <- degparms_optim_names + n_plot_rows = length(obs_vars) + n_plots = n_plot_rows * 2 - odeini_optim <- degparms_optim[odeini_optim_names] - names(odeini_optim) <- gsub("_0$", "", names(odeini_optim)) - - odeparms_optim_trans <- degparms_optim[odeparms_optim_names] - odeparms_optim <- backtransform_odeparms(odeparms_optim_trans, - fit_1$mkinmod, - transform_rates = fit_1$transform_rates, - transform_fractions = fit_1$transform_fractions) + # Set relative plot heights, so the first and the last plot are the norm + # and the middle plots (if n_plot_rows >2) are smaller by rel.height.middle + rel.heights <- if (n_plot_rows > 2) { + c(1, rep(rel.height.middle, n_plot_rows - 2), 1) + } else { + rep(1, n_plot_rows) + } - fit_a <- x$mmkin_orig[[a]] + layout_matrix = matrix(1:n_plots, + n_plot_rows, 2, byrow = TRUE) + layout(layout_matrix, heights = rel.heights) - state_ini <- fit_a$bparms.state - state_ini[names(odeini_optim)] <- odeini_optim + solution_type = fit_1$solution_type - odeparms <- fit_a$bparms.ode - odeparms[names(odeparms_optim)] <- odeparms_optim + outtimes <- sort(unique(c(x$data$time, + seq(xlim[1], xlim[2], length.out = 50)))) - mkinfit_call[["observed"]] <- ds[[a]] - mkinfit_call[["parms.ini"]] <- odeparms - mkinfit_call[["state.ini"]] <- state_ini + pred_ds <- purrr::map_dfr(i, function(ds_i) { + odeparms_trans <- degparms_all[ds_i, odeparms_names] + odeparms <- backtransform_odeparms(odeparms_trans, + x$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) - mkinfit_call[["control"]] <- list(iter.max = 0) - mkinfit_call[["quiet"]] <- TRUE + odeini <- degparms_all[ds_i, odeini_names] + names(odeini) <- gsub("_0", "", names(odeini)) - res <- suppressWarnings(do.call("mkinfit", mkinfit_call)) - return(res) + 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])) }) - # Set dimensions with names and the class (mmkin) - attributes(mmkin_nlme) <- attributes(x$mmkin_orig[, i]) - - plot(mmkin_nlme, main = main, legends = legends, - resplot = resplot, standardized = standardized, - show_errmin = show_errmin, - errmin_var = errmin_var, errmin_digits = errmin_digits, - cex = cex, - rel.height.middle = rel.height.middle, - ymax = ymax, ...) - + degparms_all_pop <- c(fixef(x), 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", "", names(odeini_pop)) + + 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 top row of plots when we have more than one row + # Reduce bottom margin by 2.1 - hides x axis legend + if (plot_row == 1 & n_plot_rows > 1) { + par(mar = c(3.0, 4.1, 4.1, 2.1)) + } + + # Margins for middle rows of plots, if any + if (plot_row > 1 & plot_row < n_plot_rows) { + # Reduce top margin by 2 after the first plot as we have no main title, + # reduced plot height, therefore we need rel.height.middle in the layout + par(mar = c(3.0, 4.1, 2.1, 2.1)) + } + + # Margins for bottom row of plots when we have more than one row + if (plot_row == n_plot_rows & n_plot_rows > 1) { + # Restore bottom margin for last plot to show x axis legend + par(mar = c(5.1, 4.1, 2.1, 2.1)) + } + + plot(pred_pop$time, pred_pop[[obs_var]], + main = obs_var, + type = "l", lwd = 2, + xlim = xlim, ylim = ylim_row, + xlab = xlab, ylab = ylab, frame = frame, + cex = cex) + + 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 (plot_row %in% legends) { + legend(lpos, inset = inset, + legend = c("Population mean", 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)) + } + + 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", + main = obs_var, + 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", + main = obs_var, + 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") + + 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 74abeeba..a73b8b09 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -37,6 +37,7 @@ reference: - "`[.mmkin`" - plot.mmkin - AIC.mmkin + - print.mmkin - title: Mixed models desc: Create and work with nonlinear mixed models contents: @@ -44,6 +45,7 @@ reference: - plot.nlme.mmkin - nlme_function - get_deg_func + - saemix_model - title: Datasets and known results contents: - FOCUS_2006_A diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html index 190bfe87..c29e413d 100644 --- a/docs/dev/news/index.html +++ b/docs/dev/news/index.html @@ -146,7 +146,9 @@ mkin 0.9.50.4 (unreleased) Unreleased
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index 4bca36e6..8c97cd7c 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-10-22T11:16Z +last_built: 2020-10-23T22:30Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/dev/reference/Extract.mmkin.html b/docs/dev/reference/Extract.mmkin.html index c69259b6..0c02355f 100644 --- a/docs/dev/reference/Extract.mmkin.html +++ b/docs/dev/reference/Extract.mmkin.html @@ -40,7 +40,7 @@ - + @@ -72,7 +72,7 @@ mkin - 0.9.50.3 + 0.9.50.4
@@ -120,7 +120,7 @@