aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2016-06-25 19:11:21 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2016-06-25 19:11:21 +0200
commite8392a8e110bb1957adc9e2047642f9387ff83db (patch)
tree013c8ccd7a821eef12caba24d395982822c8e417 /R
parent5f4a25fad9a5323611855145e6b31267b3ed9e50 (diff)
Working state, but not backwards compatible
In this commit, the separatation of plots for observed variables works. The formatting should be improved. However, in gmkin, plotting with show_residuals = TRUE is used, and the GUI expects that the residual plot is below the main plot. This behaviour should be restored.
Diffstat (limited to 'R')
-rw-r--r--R/plot.mkinfit.R136
1 files changed, 88 insertions, 48 deletions
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)
}

Contact - Imprint