Plot more stresses and effects. Closes #21 and fixes parts of #7

NEWS.md

 # stressaddition (development version)
 
-* Added predict_mixture() for the prediction of the effects of mixtures of two toxicants.
-* Fixed documentation of `ecxsys()` and `predict_ecxsys`.
+* The functions `plot_effect()` and `plot_stress()` gain a `which` argument that controls which curves are plotted. Consequently, the `show_LL5_model` argument of `plot_effect()` was removed.
+* Added arguments `xlab` and `ylab` to `plot_stress`.
+* Changed some colors of the stress curves so they better match with the colors of related effect curves.
+* Added `predict_mixture()` for the prediction of the effects of mixtures of two toxicants.
+* Fixed documentation of `ecxsys()` and `predict_ecxsys()`.
 
 # stressaddition 2.0.0
 

R/helpers.R

@@ -33,12 +33,12 @@ get_log_ticks <- function(x) {
 }
 
 
-adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
+adjust_smooth_concentrations <- function(model) {
     # Helper for the plot functions, not exported.
     # Deals with the concentrations which are unnecessary for plotting.
     # This means it removes the concentrations in the gap and scales the
     # concentrations left of the gap up.
-
+    curves <- model$curves
     gap_idx <- min(which(!curves$use_for_plotting))
 
     # Keep only the values to the left and right of the axis break:
@@ -50,7 +50,7 @@ adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
 
     # Shift the small concentrations upwards so the plot has a nice x-axis:
     curves$concentration[1:gap_idx] <-
-        curves$concentration[1:gap_idx] / conc_adjust_factor
+        curves$concentration[1:gap_idx] / model$conc_adjust_factor
 
     list(
         curves = curves,

R/plot_ecxsys.R

@@ -3,16 +3,18 @@
 
 #' Plot the results of the ECx-SyS model
 #'
-#' Convenience functions to plot the observed and modeled effect and the
-#' system stress with and without environmental stress.
+#' Plot the observed and modeled effects and stresses.
 #'
 #' @name plot_ecxsys
 #'
 #' @param model The list returned from \code{\link{ecxsys}}.
-#' @param show_LL5_model Should the log-logistic models be plotted for
-#'   comparison? Defaults to \code{FALSE}.
-#' @param show_legend Should the plot include a legend? Defaults to FALSE
-#'   because it may cover some points or lines depending on the plot size.
+#' @param which A vector of curve names to plot. Allowed values are the column
+#'   names of the \code{model$curves} data frame. The default \code{NULL} only
+#'   plots the most important curves. Use \code{which = "all"} to display all
+#'   curves.
+#' @param show_legend Should the plot include a legend? Defaults to \code{FALSE}
+#'   because it may cover some parts of the plot depending on the plot size and
+#'   the number of elements shown.
 #' @param xlab,ylab Axis labels.
 #'
 #' @examples model <- ecxsys(
@@ -23,4 +25,8 @@
 #' )
 #' plot_effect(model)
 #' plot_stress(model)
+#'
+#' # Plot all curves:
+#' plot_effect(model, which = "all")
+#' plot_stress(model, which = "all")
 NULL

R/plot_effect.R

@@ -2,33 +2,47 @@
 #' @export
 plot_effect <- function(model,
                         which = NULL,
-                        show_observed = TRUE,
                         show_legend = FALSE,
                         xlab = "concentration",
                         ylab = "effect") {
     stopifnot(inherits(model, "ecxsys"))
+
     curve_names <- names(model$curves)
     valid_names <- curve_names[startsWith(curve_names, "effect")]
     if (is.null(which)) {
-        which <- c("effect_tox_sys", "effect_tox_env_sys")
-    } else if (length(which) == 1 && which == "all") {
-        which <- valid_names
+        which <- c("effect_tox_sys")
+        if (model$with_env) {
+            which <- c(which, "effect_tox_env_sys")
+        }
+    } else if ("all" %in% which) {
+        if (length(which) == 1) {
+            which <- valid_names
+        } else {
+            stop("'all' must not be combined with other curve names.")
+        }
     } else if (!all(which %in% valid_names)) {
-        stop("Argument 'which' contains invalid names. Only effects are allowed.")
-    }
-    if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
-        warning("'which' contains names with 'env' but the model was built ",
-                "without environmental effects.")
+        warning("Argument 'which' contains invalid names.")
+        if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
+            warning("'which' contains names with 'env' but the model was",
+                    " built without environmental effects.")
+        }
+        which <- which[which %in% valid_names]
     }
 
-    temp <- adjust_smooth_concentrations(
-        model$curves,
-        model$conc_adjust_factor
-    )
+    temp <- adjust_smooth_concentrations(model)
     curves <- temp$curves
     log_ticks <- get_log_ticks(curves$concentration)
     concentration <- c(curves$concentration[1], model$args$concentration[-1])
 
+    legend_df <- data.frame(
+        text = character(),
+        pch = numeric(),
+        lty = numeric(),
+        col = character(),
+        order = numeric(),  # controls sorting of legend elements
+        stringsAsFactors = FALSE
+    )
+
     plot(
         NA,
         NA,
@@ -42,37 +56,31 @@ plot_effect <- function(model,
         las = 1
     )
 
-    if (show_observed) {
+    points(
+        concentration,
+        model$args$effect_tox_observed,
+        pch = 16,
+        col = "blue"
+    )
+    legend_df[nrow(legend_df) + 1, ] <- list("tox (observed)", 16, 0, "blue", 1)
+    if (model$with_env) {
         points(
             concentration,
-            model$args$effect_tox_observed,
+            model$args$effect_tox_env_observed,
             pch = 16,
-            col = "blue"
-        )
-        if (model$with_env) {
-            points(
-                concentration,
-                model$args$effect_tox_env_observed,
-                pch = 16,
-                col = "red"
-            )
-        }
-    }
-    if ("effect_tox_LL5" %in% which) {
-        lines(
-            curves$concentration,
-            curves$effect_tox_LL5,
-            col = "darkblue",
-            lty = 4
+            col = "red"
         )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox + env (observed)", 16, 0, "red", 2)
     }
-
+    # The lines are drawn in this order to ensure that dotted and dashed lines
+    # are on top of solid lines for better visibility.
     if ("effect_tox_sys" %in% which) {
         lines(
             curves$concentration,
             curves$effect_tox_sys,
             col = "blue"
         )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox + sys", NA, 1, "blue", 5)
     }
     if ("effect_tox" %in% which) {
         lines(
@@ -81,22 +89,25 @@ plot_effect <- function(model,
             col = "deepskyblue",
             lty = 2
         )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox", NA, 2, "deepskyblue", 4)
+    }
+    if ("effect_tox_LL5" %in% which) {
+        lines(
+            curves$concentration,
+            curves$effect_tox_LL5,
+            col = "darkblue",
+            lty = 3
+        )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox (LL5)", NA, 3, "darkblue", 3)
     }
     if (model$with_env) {
-        if ("effect_tox_env_LL5" %in% which) {
-            lines(
-                curves$concentration,
-                curves$effect_tox_env_LL5,
-                col = "darkred",
-                lty = 4
-            )
-        }
         if ("effect_tox_env_sys" %in% which) {
             lines(
                 curves$concentration,
                 curves$effect_tox_env_sys,
                 col = "red"
             )
+            legend_df[nrow(legend_df) + 1, ] <- list("tox + env + sys", NA, 1, "red", 8)
         }
         if ("effect_tox_env" %in% which) {
             lines(
@@ -105,6 +116,16 @@ plot_effect <- function(model,
                 col = "orange",
                 lty = 2
             )
+            legend_df[nrow(legend_df) + 1, ] <- list("tox + env", NA, 2, "orange", 7)
+        }
+        if ("effect_tox_env_LL5" %in% which) {
+            lines(
+                curves$concentration,
+                curves$effect_tox_env_LL5,
+                col = "darkred",
+                lty = 3
+            )
+            legend_df[nrow(legend_df) + 1, ] <- list("tox + env (LL5)", NA, 3, "darkred", 6)
         }
     }
 
@@ -113,28 +134,13 @@ plot_effect <- function(model,
     plotrix::axis.break(1, breakpos = temp$axis_break_conc)
 
     if (show_legend) {
-        legend_text <- c("tox", "tox + sys")
-        legend_pch <- c(NA, 16)
-        legend_lty <- c(2, 1)
-        legend_col <- c("deepskyblue", "blue")
-        if (model$with_env) {
-            legend_text <- c(legend_text, "tox + env", "tox + env + sys")
-            legend_pch <- c(legend_pch, NA, 16)
-            legend_lty <- c(legend_lty, 2, 1)
-            legend_col <- c(legend_col, "orange", "red")
-        }
-        if (show_LL5_model) {
-            legend_text <- c(legend_text, "tox (LL5)", "tox + env (LL5)")
-            legend_pch <- c(legend_pch, NA, NA)
-            legend_lty <- c(legend_lty, 4, 4)
-            legend_col <- c(legend_col, "darkblue", "darkred")
-        }
+        legend_df <- legend_df[order(legend_df$order), ]
         legend(
             "topright",
-            legend = legend_text,
-            pch = legend_pch,
-            lty = legend_lty,
-            col = legend_col
+            legend = legend_df$text,
+            pch = legend_df$pch,
+            lty = legend_df$lty,
+            col = legend_df$col
         )
     }
 }

R/plot_stress.R

 #' @rdname plot_ecxsys
 #' @export
-plot_stress <- function(model, show_legend = FALSE) {
+plot_stress <- function(model,
+                        which = NULL,
+                        show_legend = FALSE,
+                        xlab = "concentration",
+                        ylab = "stress") {
     stopifnot(inherits(model, "ecxsys"))
-    temp <- adjust_smooth_concentrations(
-        model$curves,
-        model$conc_adjust_factor
-    )
+
+    curve_names <- names(model$curves)
+    valid_names <- curve_names[
+        startsWith(curve_names, "stress") | startsWith(curve_names, "sys")
+    ]
+    if (is.null(which)) {
+        which <- c("sys_tox")
+        if (model$with_env) {
+            which <- c(which, "sys_tox_env")
+        }
+    } else if ("all" %in% which) {
+        if (length(which) == 1) {
+            which <- valid_names
+        } else {
+            stop("'all' must not be combined with other curve names.")
+        }
+    } else if (!all(which %in% valid_names)) {
+        warning("Argument 'which' contains invalid names.")
+        if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
+            warning("'which' contains names with 'env' but the model was",
+                    " built without environmental effects.")
+        }
+        which <- which[which %in% valid_names]
+        if (length(which) == 0) {
+            stop("No curves to display.")
+        }
+    }
+
+    temp <- adjust_smooth_concentrations(model)
     curves <- temp$curves
     axis_break_conc <- temp$axis_break_conc
     log_ticks <- get_log_ticks(curves$concentration)
     concentration <- c(curves$concentration[1], model$args$concentration[-1])
 
+    legend_df <- data.frame(
+        text = character(),
+        pch = numeric(),
+        lty = numeric(),
+        col = character(),
+        order = numeric(),  # controls sorting of legend elements
+        stringsAsFactors = FALSE
+    )
+
     plot(
         NA,
         NA,
         xlim = range(curves$concentration, na.rm = TRUE),
-        ylim = c(0, 1),
+        ylim = c(0, max(curves[, which], 1, na.rm = TRUE)),
         log = "x",
-        xlab = "concentration",
-        ylab = "system stress",
+        xlab = xlab,
+        ylab = ylab,
         xaxt = "n",
         las = 1,
         bty = "L"
     )
 
-    lines(
-        curves$concentration,
-        curves$sys_tox,
-        col = "blue"
-    )
-    points(
-        concentration,
-        model$sys_tox_not_fitted,
-        pch = 16,
-        col = "blue"
-    )
-
-    if (model$with_env) {
+    # The lines are drawn in this order to ensure that dotted and dashed lines
+    # are on top of solid lines for better visibility.
+    if ("stress_tox_sys" %in% which) {
+        lines(
+            curves$concentration,
+            curves$stress_tox_sys,
+            col = "blue"
+        )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox + sys", NA, 1, "blue", 3)
+    }
+    if ("stress_tox" %in% which) {
         lines(
             curves$concentration,
-            curves$sys_tox_env,
-            col = "red"
+            curves$stress_tox,
+            col = "deepskyblue",
+            lty = 2
         )
+        legend_df[nrow(legend_df) + 1, ] <- list("tox", NA, 2, "deepskyblue", 1)
+    }
+    if ("sys_tox" %in% which) {
         points(
             concentration,
-            model$sys_tox_env_not_fitted,
+            model$sys_tox_not_fitted,
             pch = 16,
-            col = "red"
+            col = "steelblue3"
         )
+        lines(
+            curves$concentration,
+            curves$sys_tox,
+            col = "steelblue3"
+        )
+        legend_df[nrow(legend_df) + 1, ] <- list("sys (tox)", 16, 1, "steelblue3", 2)
+    }
+    if (model$with_env) {
+        if ("stress_tox_env_sys" %in% which) {
+            lines(
+                curves$concentration,
+                curves$stress_tox_env_sys,
+                col = "red"
+            )
+            legend_df[nrow(legend_df) + 1, ] <- list("tox + env + sys", NA, 1, "red", 7)
+        }
+        if ("stress_env" %in% which) {
+            lines(
+                curves$concentration,
+                curves$stress_env,
+                col = "forestgreen",
+                lty = 3
+            )
+            legend_df[nrow(legend_df) + 1, ] <- list("env", NA, 3, "forestgreen", 4)
+        }
+        if ("stress_tox_env" %in% which) {
+            lines(
+                curves$concentration,
+                curves$stress_tox_env,
+                col = "orange",
+                lty = 2
+            )
+            legend_df[nrow(legend_df) + 1, ] <- list("tox + env", NA, 2, "orange", 5)
+        }
+        if ("sys_tox_env" %in% which) {
+            points(
+                concentration,
+                model$sys_tox_env_not_fitted,
+                pch = 16,
+                col = "violetred"
+            )
+            lines(
+                curves$concentration,
+                curves$sys_tox_env,
+                col = "violetred"
+            )
+            legend_df[nrow(legend_df) + 1, ] <- list("sys (tox + env)", 16, 1, "violetred", 6)
+        }
     }
 
     axis(1, at = log_ticks$major, labels = log_ticks$major_labels)
@@ -55,18 +143,13 @@ plot_stress <- function(model, show_legend = FALSE) {
     plotrix::axis.break(1, breakpos = axis_break_conc)
 
     if (show_legend) {
-        legend_text <- c("tox")
-        legend_col <- c("blue")
-        if (model$with_env) {
-            legend_text <- c(legend_text, "tox + env")
-            legend_col <- c(legend_col, "red")
-        }
+        legend_df <- legend_df[order(legend_df$order), ]
         legend(
             "topright",
-            legend = legend_text,
-            col = legend_col,
-            pch = 16,
-            lty = 1
+            legend = legend_df$text,
+            pch = legend_df$pch,
+            lty = legend_df$lty,
+            col = legend_df$col
         )
     }
 }

man/plot_ecxsys.Rd

@@ -10,7 +10,7 @@ plot_effect(
   model,
   which = NULL,
   show_observed = TRUE,
-  show_legend = FALSE,
+  show_legend = TRUE,
   xlab = "concentration",
   ylab = "effect"
 )
@@ -20,17 +20,20 @@ plot_stress(model, show_legend = FALSE)
 \arguments{
 \item{model}{The list returned from \code{\link{ecxsys}}.}
 
-\item{show_legend}{Should the plot include a legend? Defaults to FALSE
-because it may cover some points or lines depending on the plot size.}
+\item{which}{A vector of curve names to plot. Allowed values are the column
+names of the \code{model$curves} data frame. The default \code{NULL} only
+plots the most important curves. Use \code{which = "all"} to display all
+curves.}
 
-\item{xlab, ylab}{Axis labels.}
+\item{show_observed}{Should the observed values be plotted as points?}
+
+\item{show_legend}{Should the plot include a legend? May cover some parts of
+the plot depending on the plot size and the number of elements shown.}
 
-\item{show_LL5_model}{Should the log-logistic models be plotted for
-comparison? Defaults to \code{FALSE}.}
+\item{xlab, ylab}{Axis labels.}
 }
 \description{
-Convenience functions to plot the observed and modeled effect and the
-system stress with and without environmental stress.
+Plot the observed and modeled effects and stresses.
 }
 \examples{
 model <- ecxsys(
@@ -41,4 +44,8 @@ model <- ecxsys(
 )
 plot_effect(model)
 plot_stress(model)
+
+# Plot all curves:
+plot_effect(model, which = "all")
+plot_stress(model, which = "all")
 }