From a1631098acfc3352e19c331e568bd6f5766b3c3d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 9 Nov 2020 16:33:19 +0100 Subject: Merge plot methods for nlme.mmkin and saem.mmkin This avoids code duplication --- NAMESPACE | 2 +- R/nlme.mmkin.R | 3 + R/plot.nlme.mmkin.R | 216 --------------------- R/plot.saem.mmkin.R | 218 --------------------- R/plot_mixed.R | 294 +++++++++++++++++++++++++++++ docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Rplot002.png | Bin 56909 -> 57363 bytes docs/dev/reference/Rplot003.png | Bin 52903 -> 56909 bytes docs/dev/reference/dot-plot_mixed.html | 203 ++++++++++++++++++++ docs/dev/reference/get_deg_func.html | 10 +- docs/dev/reference/index.html | 10 +- docs/dev/reference/plot_mixed-1.png | Bin 0 -> 86076 bytes docs/dev/reference/plot_mixed-2.png | Bin 0 -> 164488 bytes docs/dev/reference/plot_mixed-3.png | Bin 0 -> 164014 bytes docs/dev/reference/plot_mixed.html | 336 +++++++++++++++++++++++++++++++++ docs/dev/reference/reexports.html | 14 +- docs/dev/sitemap.xml | 5 +- man/plot.nlme.mmkin.Rd | 94 --------- man/plot.saem.mmkin.Rd | 94 --------- man/plot_mixed.Rd | 123 ++++++++++++ man/reexports.Rd | 5 +- 21 files changed, 982 insertions(+), 647 deletions(-) delete mode 100644 R/plot.nlme.mmkin.R delete mode 100644 R/plot.saem.mmkin.R create mode 100644 R/plot_mixed.R create mode 100644 docs/dev/reference/dot-plot_mixed.html create mode 100644 docs/dev/reference/plot_mixed-1.png create mode 100644 docs/dev/reference/plot_mixed-2.png create mode 100644 docs/dev/reference/plot_mixed-3.png create mode 100644 docs/dev/reference/plot_mixed.html delete mode 100644 man/plot.nlme.mmkin.Rd delete mode 100644 man/plot.saem.mmkin.Rd create mode 100644 man/plot_mixed.Rd diff --git a/NAMESPACE b/NAMESPACE index 29462082..e6cbe2a7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -74,6 +74,7 @@ export(mkinresplot) export(mkinsub) export(mmkin) export(nafta) +export(nlme) export(nlme_data) export(nlme_function) export(parms) @@ -99,7 +100,6 @@ 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/nlme.mmkin.R b/R/nlme.mmkin.R index ca849b85..695c63e9 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -6,6 +6,9 @@ # following environment will be in the mkin namespace. .nlme_env <- new.env(parent = emptyenv()) +#' @export +nlme::nlme + #' Retrieve a degradation function from the mmkin namespace #' #' @importFrom utils getFromNamespace diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R deleted file mode 100644 index 4228b18a..00000000 --- a/R/plot.nlme.mmkin.R +++ /dev/null @@ -1,216 +0,0 @@ -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 -#' @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 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) 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) -#' -#' library(nlme) -#' # 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) -#' } -#' @export -plot.nlme.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 <- 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), - 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, "pearson", "response") - - observed <- cbind(x$data, - residual = residuals(x, type = 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_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", "", 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/R/plot.saem.mmkin.R b/R/plot.saem.mmkin.R deleted file mode 100644 index ce43fdb6..00000000 --- a/R/plot.saem.mmkin.R +++ /dev/null @@ -1,218 +0,0 @@ -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/R/plot_mixed.R b/R/plot_mixed.R new file mode 100644 index 00000000..68404de4 --- /dev/null +++ b/R/plot_mixed.R @@ -0,0 +1,294 @@ +if(getRversion() >= '2.15.1') utils::globalVariables("ds") + +#' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object +#' +#' @name plot_mixed +#' @param x An object of class [saem.mmkin] or [nlme.mmkin] +#' @param i A numeric index to select datasets for which to plot the individual predictions, +#' 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 stats coefficients +#' @return The functions are called for their 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) +#' +#' library(nlme) +#' # 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) +#' plot(f_saem) +#' } +#' @rdname plot_mixed +#' @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, ...) +{ + fit_1 <- x$mmkin[[1]] + ds_names <- colnames(x$mmkin) + + degparms_optim <- saemix::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 + + residual_type = ifelse(standardized, "iwres", "ires") + + residuals <- x$so@results@predictions[[residual_type]] + + degparms_pop <- x$so@results@fixed.effects + names(degparms_pop) <- degparms_optim_names + + .plot_mixed(x, i, + degparms_optim, degparms_pop, residuals, + obs_vars, standardized, xlab, xlim, + resplot, ymax, maxabs, + ncol.legend, nrow.legend, + rel.height.legend, rel.height.bottom, + pch_ds, col_ds, lty_ds, frame, ...) +} + +#' @rdname plot_mixed +#' @export +plot.nlme.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, ...) +{ + degparms_optim <- coefficients(x) + degparms_pop <- nlme::fixef(x) + + residuals <- residuals(x, + type = ifelse(standardized, "pearson", "response")) + + .plot_mixed(x, i, + degparms_optim, degparms_pop, residuals, + obs_vars, standardized, xlab, xlim, + resplot, ymax, maxabs, + ncol.legend, nrow.legend, + rel.height.legend, rel.height.bottom, + pch_ds, col_ds, lty_ds, frame, ...) +} + +.plot_mixed <- function(x, i = 1:ncol(x$mmkin), + degparms_optim, + degparms_pop, + residuals, + 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_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(names(degparms_optim), names(degparms_fixed)) + colnames(degparms_all) <- degparms_all_names + + degparms_all_pop <- c(degparms_pop, degparms_fixed) + + 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 = residuals) + + 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])) + }) + + 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/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index e14af5e0..4b2b76f7 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-09T13:23Z +last_built: 2020-11-09T15:29Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/dev/reference/Rplot002.png b/docs/dev/reference/Rplot002.png index cd2014eb..8ada7133 100644 Binary files a/docs/dev/reference/Rplot002.png and b/docs/dev/reference/Rplot002.png differ diff --git a/docs/dev/reference/Rplot003.png b/docs/dev/reference/Rplot003.png index 7514f1a4..cd2014eb 100644 Binary files a/docs/dev/reference/Rplot003.png and b/docs/dev/reference/Rplot003.png differ diff --git a/docs/dev/reference/dot-plot_mixed.html b/docs/dev/reference/dot-plot_mixed.html new file mode 100644 index 00000000..e359c5ea --- /dev/null +++ b/docs/dev/reference/dot-plot_mixed.html @@ -0,0 +1,203 @@ + + + + + + + + +Not intended to be called directly by the user — .plot_mixed • mkin + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Not intended to be called directly by the user

+ Source: R/plot_mixed.R + +
+ +
+

Not intended to be called directly by the user

+
+ + 
.plot_mixed(
+  x,
+  i = 1:ncol(x$mmkin),
+  degparms_optim,
+  degparms_pop,
+  residuals,
+  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,
+  ...
+)
+ + + +
+ +
+ + +
+
+

Developed by Johannes Ranke.

+
+ +
+

Site built with pkgdown 1.6.1.

+
+ +
+
+ + + + + + + + diff --git a/docs/dev/reference/get_deg_func.html b/docs/dev/reference/get_deg_func.html index ea0676cc..a266bf5f 100644 --- a/docs/dev/reference/get_deg_func.html +++ b/docs/dev/reference/get_deg_func.html @@ -72,7 +72,7 @@ mkin - 0.9.50.3 + 0.9.50.4 @@ -120,7 +120,7 @@