diff options
| -rw-r--r-- | NEWS.md | 8 | ||||
| -rw-r--r-- | R/plot.mkinfit.R | 136 | ||||
| -rw-r--r-- | man/mmkin.Rd | 8 | ||||
| -rw-r--r-- | man/plot.mkinfit.Rd | 49 | ||||
| -rw-r--r-- | man/plot.mmkin.Rd | 2 |
5 files changed, 134 insertions, 69 deletions
@@ -6,6 +6,10 @@ - The title was changed to `Kinetic evaluations of chemical degradation data` +- `plot.mkinfit`: If a residual plot is requested, show it next to the plot of the fit, not below. This may cause existing code to produce bad-looking plots, but was done to make the next feature possible without increasing code complexity too much. + +- `plot.mkinfit`: Add the possibility to show fits and residual plots separately for the observed variables + - The main vignette `mkin` was converted to R markdown and updated - The function `add_err` was added to the package, making it easy to generate simulated data using an error model based on the normal distribution @@ -26,7 +30,9 @@ - `mkinfit`: Check for all observed variables when checking if the user tried to fix formation fractions when fitting them using ilr transformation. -- `plot.mmkin`: Removed some leftover code that set the plot margins wrongly in the case of a single fit to be plotted, so the main title was misplaced +- `plot.mmkin`: Removed some leftover code that set the plot margins wrongly in the case of a single fit to be plotted, so the main title was misplaced. + +- `plot.mkinfit`: Correct default values for `col_obs`, `pch_obs` and `lty_obs` for the case that `obs_vars` is specified. ## mkin 0.9.42 (2016-03-25) diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 6ca7e531..6bc64351 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2015 Johannes Ranke +# Copyright (C) 2010-2016 Johannes Ranke # Contact: jranke@uni-bremen.de # This file is part of the R package mkin @@ -22,18 +22,16 @@ plot.mkinfit <- function(x, fit = x, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = "default", - col_obs = 1:length(fit$mkinmod$map), + col_obs = 1:length(obs_vars), pch_obs = col_obs, - lty_obs = rep(1, length(fit$mkinmod$map)), + lty_obs = rep(1, length(obs_vars)), add = FALSE, legend = !add, show_residuals = FALSE, maxabs = "auto", + sep_obs = FALSE, rel.height.middle = 0.9, lpos = "topright", inset = c(0.05, 0.05), ...) { if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") - - if (ylim[[1]] == "default") { - ylim = c(0, max(subset(fit$data, variable %in% obs_vars)$observed, na.rm = TRUE)) - } + if (add && sep_obs) stop("If adding to an existing plot we can not show observed variables separately") solution_type = fit$solution_type parms.all <- c(fit$bparms.optim, fit$bparms.fixed) @@ -57,48 +55,90 @@ plot.mkinfit <- function(x, fit = x, out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol) - # Set up the plot if not to be added to an existing plot - if (add == FALSE) { - if (show_residuals) { - oldpar <- par(no.readonly = TRUE) - layout(matrix(c(1, 2), 2, 1), heights = c(2, 1.3)) - par(mar = c(3, 4, 4, 2) + 0.1) - } - plot(0, type="n", - xlim = xlim, ylim = ylim, - xlab = xlab, ylab = ylab, ...) - } - # Plot the data and model output - names(col_obs) <- names(pch_obs) <- names(lty_obs) <- names(fit$mkinmod$map) - for (obs_var in obs_vars) { - points(subset(fit$data, variable == obs_var, c(time, observed)), - pch = pch_obs[obs_var], col = col_obs[obs_var]) - } - matlines(out$time, out[obs_vars], col = col_obs[obs_vars], lty = lty_obs[obs_vars]) - if (legend == TRUE) { - legend_names = lapply(names(fit$mkinmod$spec), function(x) { - if (!is.null(fit$mkinmod$spec[[x]]$full_name)) - if (is.na(fit$mkinmod$spec[[x]]$full_name)) x - else fit$mkinmod$spec[[x]]$full_name - else x - }) - legend(lpos, inset= inset, legend = legend_names, - col = col_obs[obs_vars], pch = pch_obs[obs_vars], lty = lty_obs[obs_vars]) + names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars + + # Create a plot layout only if not to be added to an existing plot + # or only a single plot is requested (e.g. by plot.mmkin) + do_layout = FALSE + if (show_residuals) do_layout = TRUE + if (sep_obs) do_layout = TRUE + n_plot_rows = if (sep_obs) length(obs_vars) else 1 + + if (do_layout) { + # Layout should be restored afterwards + oldpar <- par(no.readonly = TRUE) + + n_plot_cols = if (show_residuals) 2 else 1 + n_plots = n_plot_rows * n_plot_cols + + # 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) + layout_matrix = matrix(1:n_plots, + n_plot_rows, n_plot_cols, byrow = TRUE) + layout(layout_matrix, heights = rel.heights) } - # Show residuals if requested - if (show_residuals) { - par(mar = c(5, 4, 0, 2) + 0.1) - residuals <- subset(fit$data, variable %in% obs_vars, residual) - if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) - plot(0, type="n", - xlim = xlim, - ylim = c(-1.2 * maxabs, 1.2 * maxabs), - xlab = xlab, ylab = "Residuals") - for(obs_var in obs_vars){ - residuals_plot <- subset(fit$data, variable == obs_var, c("time", "residual")) - points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var]) + + # Replicate legend position argument if necessary + if (length(lpos) == 1) lpos = rep(lpos, n_plot_rows) + + # Loop over plot rows + for (plot_row in 1:n_plot_rows) { + + row_obs_vars = if (sep_obs) obs_vars[plot_row] else obs_vars + + # Set ylim to sensible default, or to the specified value + if (ylim[[1]] == "default") { + ylim_row = c(0, max(c(subset(fit$data, variable %in% row_obs_vars)$observed, + subset(fit$data, variable %in% row_obs_vars)$fitted), + na.rm = TRUE)) + } else { + ylim_row = ylim + } + + # Set up the main plot if not to be added to an existing plot + if (add == FALSE) { + plot(0, type="n", + xlim = xlim, ylim = ylim_row, + xlab = xlab, ylab = ylab, ...) + } + + # Plot the data + for (obs_var in row_obs_vars) { + points(subset(fit$data, variable == obs_var, c(time, observed)), + pch = pch_obs[obs_var], col = col_obs[obs_var]) + } + + # Plot the model output + matlines(out$time, out[row_obs_vars], col = col_obs[row_obs_vars], lty = lty_obs[row_obs_vars]) + + if (legend == TRUE) { + # Get full names from model definition if they are available + legend_names = lapply(row_obs_vars, function(x) { + if (!is.null(fit$mkinmod$spec[[x]]$full_name)) + if (is.na(fit$mkinmod$spec[[x]]$full_name)) x + else fit$mkinmod$spec[[x]]$full_name + else x + }) + legend(lpos[plot_row], inset= inset, legend = legend_names, + col = col_obs[row_obs_vars], pch = pch_obs[row_obs_vars], lty = lty_obs[row_obs_vars]) + } + + # Show residuals if requested + if (show_residuals) { + residuals <- subset(fit$data, variable %in% row_obs_vars, residual) + if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) + plot(0, type="n", + xlim = xlim, + ylim = c(-1.2 * maxabs, 1.2 * maxabs), + xlab = xlab, ylab = "Residuals") + for(obs_var in row_obs_vars){ + residuals_plot <- subset(fit$data, variable == obs_var, c("time", "residual")) + points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var]) + } + abline(h = 0, lty = 2) } - abline(h = 0, lty = 2) - par(oldpar, no.readonly = TRUE) } + if (do_layout) par(oldpar, no.readonly = TRUE) } diff --git a/man/mmkin.Rd b/man/mmkin.Rd index f701890a..44caf8d2 100644 --- a/man/mmkin.Rd +++ b/man/mmkin.Rd @@ -56,13 +56,19 @@ m_synth_FOMC_lin <- mkinmod(parent = list(type = "FOMC", to = "M1"), models <- list(SFO_lin = m_synth_SFO_lin, FOMC_lin = m_synth_FOMC_lin) datasets <- lapply(synthetic_data_for_UBA_2014[1:3], function(x) x$data) -time_default <- system.time(fits <- mmkin(models, datasets)) +time_default <- system.time(fits.0 <- mmkin(models, datasets)) time_1 <- system.time(fits.1 <- mmkin(models, datasets, cores = 1)) time_default time_1 endpoints(fits[["SFO_lin", 2]]) + +# Plot.mkinfit handles rows or columns of mmkin result objects +plot(fits.0[1, ]) +# Double brackets to select a single mkinfit object, which will be +# plotted by plot.mkinfit +plot(fits.0[[1, 1]], sep_obs = TRUE, show_residuals = TRUE) } } \keyword{ optimize } diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd index 494dc38d..b80928f7 100644 --- a/man/plot.mkinfit.Rd +++ b/man/plot.mkinfit.Rd @@ -14,10 +14,11 @@ xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = "default", - col_obs = 1:length(fit$mkinmod$map), pch_obs = col_obs, - lty_obs = rep(1, length(fit$mkinmod$map)), + col_obs = 1:length(obs_vars), pch_obs = col_obs, + lty_obs = rep(1, length(obs_vars)), add = FALSE, legend = !add, show_residuals = FALSE, maxabs = "auto", + sep_vars = FALSE, rel.height.middle = 0.9, lpos = "topright", inset = c(0.05, 0.05), \dots) } \arguments{ @@ -25,7 +26,7 @@ Alias for fit introduced for compatibility with the generic S3 method. } \item{fit}{ - an object of class \code{\link{mkinfit}}. + An object of class \code{\link{mkinfit}}. } \item{obs_vars}{ A character vector of names of the observed variables for which the @@ -33,47 +34,54 @@ in the model. } \item{xlab}{ - label for the x axis. + Label for the x axis. } \item{ylab}{ - label for the y axis. + Label for the y axis. } \item{xlim}{ - plot range in x direction. + Plot range in x direction. } \item{ylim}{ - plot range in y direction. + Plot range in y direction. } \item{col_obs}{ - colors used for plotting the observed data and the corresponding model prediction lines. + Colors used for plotting the observed data and the corresponding model prediction lines. } \item{pch_obs}{ - symbols to be used for plotting the data. + Symbols to be used for plotting the data. } \item{lty_obs}{ - line types to be used for the model predictions. + Line types to be used for the model predictions. } \item{add}{ - should the plot be added to an existing plot? + Should the plot be added to an existing plot? } \item{legend}{ - should a legend be added to the plot? + Should a legend be added to the plot? } \item{show_residuals}{ - should residuals be shown in the lower third of the plot? + Should residuals be shown in the lower third of the plot? } \item{maxabs}{ Maximum absolute value of the residuals. This is used for the scaling of the y axis and defaults to "auto". } + \item{sep_obs}{ + Should the observed variables be shown in separate subplots? + } + \item{rel.height.middle}{ + The relative height of the middle plot, if more than two rows of plots are shown. + } \item{lpos}{ - position of the legend. Passed to \code{\link{legend}} as the first argument. + Position(s) of the legend(s). Passed to \code{\link{legend}} as the first argument. + If not length one, this should be of the same length as the obs_var argument. } \item{inset}{ Passed to \code{\link{legend}} if applicable. } \item{\dots}{ - further arguments passed to \code{\link{plot}}. + Further arguments passed to \code{\link{plot}}. } } \value{ @@ -81,13 +89,18 @@ } \examples{ # One parent compound, one metabolite, both single first order, path from -# parent to sink included +# parent to sink included, use Levenberg-Marquardt for speed SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"), m1 = mkinsub("SFO", full = "Metabolite M1" )) -fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) +fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, method.modFit = "Marq") plot(fit) + +# Show the observed variables separately +plot(fit, sep_obs = TRUE) + +# Show the observed variables separately, with residuals +plot(fit, sep_obs = TRUE, show_residuals = TRUE, lpos = c("topright", "bottomright")) } \author{ Johannes Ranke } -\keyword{ hplot } diff --git a/man/plot.mmkin.Rd b/man/plot.mmkin.Rd index 8362f16c..ffd83d2f 100644 --- a/man/plot.mmkin.Rd +++ b/man/plot.mmkin.Rd @@ -32,7 +32,7 @@ Passed to the plot functions and \code{\link{mtext}}. } \item{rel.height.middle}{ - The relative height of the middle plot. + The relative height of the middle plot, if more than two rows of plots are shown. } \item{\dots}{ Further arguments passed to \code{\link{plot.mkinfit}} and \code{\link{mkinresplot}}. |