diff --git a/DESCRIPTION b/DESCRIPTION
index 422ecaa2ca550bd37402085b183fab19ea0aa1be..b396809dd5850402898b79c47bdaf32ca5c4391a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: stressaddition
 Type: Package
 Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 2.3.0
-Date: 2020-03-11
+Version: 2.4.0
+Date: 2020-03-13
 Authors@R: c(person("Sebastian", 
                     "Henz", 
                     role = c("aut", "cre"), 
@@ -31,4 +31,4 @@ Imports:
     plotrix
 Suggests: 
     testthat (>= 2.1.0)
-RoxygenNote: 7.0.2
+RoxygenNote: 7.1.0
diff --git a/NEWS.md b/NEWS.md
index 2b8b99f441931ee20633ddb96ff63ca1c671240e..f43297f4ed42b587536adf2677196cefe1837555 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,7 +1,12 @@
+# stressaddition 2.4.0
+
+* Improved `plot_effect()` and `plot_stress()`. You can now control whether the observed values (the points) should be plotted using the `which` argument.
+* Renamed `sys_tox_not_fitted` and `sys_tox_env_not_fitted` to `sys_tox_observed` and `sys_tox_env_observed` in the output of `ecxsys()`.
+
 # stressaddition 2.3.0
 
 * `predict_mixture()` accepts multiple values for the concentration of the second toxicant. Both concentration vectors must be the same length.
-* `predict_mixture()` returns a data frame with the concentrations and effects. Previously it was only a vector of effects.
+* `predict_mixture()` now returns a data frame with the concentrations and effects. Previously it was only a vector of effects.
 * `predict_mixture()` received a new argument "effect_max" which scales the returned effect values.
 * Renamed the arguments of `predict_mixture()` to use underscore letters a and b instad of 1 and 2. For example model_1 is now model_a.
 
diff --git a/R/ecxsys.R b/R/ecxsys.R
index d6ba1a7a4618efa91a3c32ed96bf5587487bda30..80e05b4f76b29e76b371c6c67dc628ba793b159c 100644
--- a/R/ecxsys.R
+++ b/R/ecxsys.R
@@ -234,7 +234,7 @@ ecxsys <- function(concentration,
         hormesis_index,
         original_options
     )
-    output$sys_tox_not_fitted <- fit_sys_output$sys[keep]
+    output$sys_tox_observed <- fit_sys_output$sys[keep]
     output$sys_tox_mod <- fit_sys_output$sys_mod
     if (inherits(fit_sys_output$sys_mod, "lm")) {
         warning(
@@ -273,7 +273,7 @@ ecxsys <- function(concentration,
             hormesis_index,
             original_options
         )
-        output$sys_tox_env_not_fitted <- fit_sys_output$sys[keep]
+        output$sys_tox_env_observed <- fit_sys_output$sys[keep]
         output$sys_tox_env_mod <- fit_sys_output$sys_mod
         if (inherits(fit_sys_output$sys_mod, "lm")) {
             warning(
diff --git a/R/plot_ecxsys.R b/R/plot_ecxsys.R
index 52a812ca66b4c61da63c89f5af78837130691d05..f1cd7ca72bf0006851b40ea7b726c556ed44e4ee 100644
--- a/R/plot_ecxsys.R
+++ b/R/plot_ecxsys.R
@@ -27,10 +27,13 @@
 #' @name plot_ecxsys
 #'
 #' @param model The list returned from \code{\link{ecxsys}}.
-#' @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 which A vector of names to plot. Allowed are the column names of the
+#'   \code{model$curves} data frame. There is also \code{"effect_tox_observed"}
+#'   and \code{"effect_tox_env_observed"} for the observed effects and
+#'   \code{"sys_tox_observed"} and \code{"sys_tox_env_observed"} for the
+#'   observed Sys. 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.
@@ -42,14 +45,18 @@
 #'     effect_tox_observed = c(90, 81, 92, 28, 0),
 #'     effect_tox_env_observed = c(29, 27, 33, 5, 0)
 #' )
-#' plot_effect(model)
-#' plot_stress(model)
+#' plot_effect(model, show_legend = TRUE)
+#' plot_stress(model, show_legend = TRUE)
 #'
 #' # Plot all curves:
 #' plot_effect(model, which = "all")
 #' plot_stress(model, which = "all")
 #'
 #' # Plot only some selected curves:
-#' plot_effect(model, which = c("effect_tox", "effect_tox_env"))
-#' plot_stress(model, which = c("stress_tox", "stress_tox_env"))
+#' plot_effect(model, which = c("effect_tox_sys", "effect_tox_env_sys"))
+#' plot_stress(model, which = c("sys_tox", "sys_tox_env"))
+#'
+#' # Plot only the observed values:
+#' plot_effect(model, which = c("effect_tox_observed", "effect_tox_env_observed"))
+#' plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed"))
 NULL
diff --git a/R/plot_effect.R b/R/plot_effect.R
index 5a0da15f0b2de26ec4e4f6bd12fd31b8b5f8542b..0645f032b95d573703537f86a828064c80e33c2e 100644
--- a/R/plot_effect.R
+++ b/R/plot_effect.R
@@ -28,19 +28,22 @@ plot_effect <- function(model,
     stopifnot(inherits(model, "ecxsys"))
 
     curve_names <- names(model$curves)
-    valid_names <- curve_names[startsWith(curve_names, "effect")]
+    valid_names <- c(
+        curve_names[startsWith(curve_names, "effect")],
+        "effect_tox_observed", "effect_tox_env_observed"  # the observed points
+    )
     if (is.null(which)) {
-        which <- c("effect_tox", "effect_tox_sys")
+        which <- c("effect_tox", "effect_tox_sys", "effect_tox_observed")
         if (model$with_env) {
-            which <- c(which, "effect_tox_env", "effect_tox_env_sys")
+            which <- c(which, "effect_tox_env", "effect_tox_env_sys",
+                       "effect_tox_env_observed")
         }
     } else if ("all" %in% which) {
-        if (length(which) == 1) {
-            which <- valid_names
-        } else {
+        if (length(which) > 1) {
             stop("'all' must not be combined with other curve names.")
         }
-    } else if (!all(which %in% valid_names)) {
+        which <- valid_names
+    } else if (any(!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",
@@ -54,15 +57,6 @@ plot_effect <- function(model,
     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,
@@ -77,31 +71,22 @@ plot_effect <- function(model,
         bty = "L"
     )
 
-    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) {
+    # 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_observed" %in% which) {
         points(
             concentration,
-            model$args$effect_tox_env_observed,
+            model$args$effect_tox_observed,
             pch = 16,
-            col = "red"
+            col = "blue"
         )
-        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(
@@ -110,7 +95,6 @@ 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(
@@ -119,16 +103,22 @@ plot_effect <- function(model,
             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_observed" %in% which) {
+            points(
+                concentration,
+                model$args$effect_tox_env_observed,
+                pch = 16,
+                col = "red"
+            )
+        }
         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(
@@ -137,7 +127,6 @@ 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(
@@ -146,7 +135,6 @@ plot_effect <- function(model,
                 col = "darkred",
                 lty = 3
             )
-            legend_df[nrow(legend_df) + 1, ] <- list("tox + env (LL5)", NA, 3, "darkred", 6)
         }
     }
 
@@ -161,13 +149,28 @@ plot_effect <- function(model,
     axis(2, col = NA, col.ticks = par("fg"), las = 1)
 
     if (show_legend) {
-        legend_df <- legend_df[order(legend_df$order), ]
-        legend(
-            "topright",
-            legend = legend_df$text,
-            pch = legend_df$pch,
-            lty = legend_df$lty,
-            col = legend_df$col
+        legend_df <- data.frame(
+            name = c("effect_tox_observed", "effect_tox", "effect_tox_sys",
+                     "effect_tox_LL5", "effect_tox_env_observed",
+                     "effect_tox_env", "effect_tox_env_sys", "effect_tox_env_LL5"),
+            text = c("tox (observed)", "tox", "tox + sys", "tox (LL5)",
+                     "tox + env (observed)", "tox + env", "tox + env + sys",
+                     "tox + env (LL5)"),
+            pch = c(16, NA, NA, NA, 16, NA, NA, NA),
+            lty = c(0, 2, 1, 3, 0, 2, 1, 3),
+            col = c("blue", "deepskyblue", "blue", "darkblue", "red", "orange",
+                    "red", "darkred"),
+            stringsAsFactors = FALSE
         )
+        legend_df <- legend_df[legend_df$name %in% which, ]
+        if (nrow(legend_df) > 0) {
+            legend(
+                "topright",
+                legend = legend_df$text,
+                pch = legend_df$pch,
+                lty = legend_df$lty,
+                col = legend_df$col
+            )
+        }
     }
 }
diff --git a/R/plot_stress.R b/R/plot_stress.R
index ee80cd49846519fae905291dd778af08b160c2c8..296f90a6b181cefe217f02aa56aebacacf4f4751 100644
--- a/R/plot_stress.R
+++ b/R/plot_stress.R
@@ -28,30 +28,27 @@ plot_stress <- function(model,
     stopifnot(inherits(model, "ecxsys"))
 
     curve_names <- names(model$curves)
-    valid_names <- curve_names[
-        startsWith(curve_names, "stress") | startsWith(curve_names, "sys")
-    ]
+    valid_names <- c(
+        curve_names[startsWith(curve_names, "stress") | startsWith(curve_names, "sys")],
+        "sys_tox_observed", "sys_tox_env_observed"  # the observed points
+    )
     if (is.null(which)) {
-        which <- c("sys_tox")
+        which <- c("sys_tox", "sys_tox_observed")
         if (model$with_env) {
-            which <- c(which, "sys_tox_env")
+            which <- c(which, "sys_tox_env", "sys_tox_env_observed")
         }
     } else if ("all" %in% which) {
-        if (length(which) == 1) {
-            which <- valid_names
-        } else {
+        if (length(which) > 1) {
             stop("'all' must not be combined with other curve names.")
         }
-    } else if (!all(which %in% valid_names)) {
+        which <- valid_names
+    } else if (any(!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_plot_concentrations(model)
@@ -59,20 +56,14 @@ plot_stress <- function(model,
     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
-    )
+    curves_w <- curves[, which[!endsWith(which, "observed")]]
+    ymax <- if (NCOL(curves_w) == 0) 1 else max(curves_w, 1, na.rm = TRUE)
 
     plot(
         NA,
         NA,
         xlim = range(curves$concentration, na.rm = TRUE),
-        ylim = c(0, max(curves[, which], 1, na.rm = TRUE)),
+        ylim = c(0, ymax),
         log = "x",
         xlab = xlab,
         ylab = ylab,
@@ -84,13 +75,20 @@ plot_stress <- function(model,
 
     # The lines are drawn in this order to ensure that dotted and dashed lines
     # are on top of solid lines for better visibility.
+    if ("sys_tox_observed" %in% which) {
+        points(
+            concentration,
+            model$sys_tox_observed,
+            pch = 16,
+            col = "steelblue3"
+        )
+    }
     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(
@@ -99,30 +97,29 @@ plot_stress <- function(model,
             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_not_fitted,
-            pch = 16,
-            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 ("sys_tox_env_observed" %in% which) {
+            points(
+                concentration,
+                model$sys_tox_env_observed,
+                pch = 16,
+                col = "violetred"
+            )
+        }
         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(
@@ -131,7 +128,6 @@ plot_stress <- function(model,
                 col = "forestgreen",
                 lty = 3
             )
-            legend_df[nrow(legend_df) + 1, ] <- list("env", NA, 3, "forestgreen", 4)
         }
         if ("stress_tox_env" %in% which) {
             lines(
@@ -140,21 +136,13 @@ plot_stress <- function(model,
                 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)
         }
     }
 
@@ -169,13 +157,28 @@ plot_stress <- function(model,
     axis(2, col = NA, col.ticks = par("fg"), las = 1)
 
     if (show_legend) {
-        legend_df <- legend_df[order(legend_df$order), ]
-        legend(
-            "topright",
-            legend = legend_df$text,
-            pch = legend_df$pch,
-            lty = legend_df$lty,
-            col = legend_df$col
+        legend_df <- data.frame(
+            name = c( "sys_tox_observed", "stress_tox", "sys_tox",
+                      "stress_tox_sys", "sys_tox_env_observed", "stress_env",
+                      "stress_tox_env", "sys_tox_env", "stress_tox_env_sys"),
+            text = c("sys (tox, observed)", "tox", "sys (tox)", "tox + sys",
+                     "sys (tox + env, observed)", "env", "tox + env",
+                     "sys (tox + env)", "tox + env + sys"),
+            pch = c(16, NA, NA, NA, 16, NA, NA, NA, NA),
+            lty = c(0, 2, 1, 1, 0, 3, 2, 1, 1),
+            col = c("steelblue3", "deepskyblue", "steelblue3", "blue",
+                    "violetred", "forestgreen", "orange", "violetred", "red"),
+            stringsAsFactors = FALSE
         )
+        legend_df <- legend_df[legend_df$name %in% which, ]
+        if (nrow(legend_df) > 0) {
+            legend(
+                "topright",
+                legend = legend_df$text,
+                pch = legend_df$pch,
+                lty = legend_df$lty,
+                col = legend_df$col
+            )
+        }
     }
 }
diff --git a/man/plot_ecxsys.Rd b/man/plot_ecxsys.Rd
index 19d5c61ed5b184d8a9d0dd157f877806c3c8b0b5..1c9c32f9517e753c541357bc979485436dd6648a 100644
--- a/man/plot_ecxsys.Rd
+++ b/man/plot_ecxsys.Rd
@@ -27,10 +27,12 @@ plot_stress(
 \arguments{
 \item{model}{The list returned from \code{\link{ecxsys}}.}
 
-\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{which}{A vector of names to plot. Allowed are the column names of the
+\code{model$curves} data frame. There is also \code{"effect_tox_observed"}
+and \code{"effect_tox_env_observed"} for the observed effects and
+\code{"sys_tox_observed"} and \code{"sys_tox_env_observed"} for the
+observed Sys. The default \code{NULL} only plots the most important curves.
+Use \code{which = "all"} to display all curves.}
 
 \item{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
@@ -48,14 +50,18 @@ model <- ecxsys(
     effect_tox_observed = c(90, 81, 92, 28, 0),
     effect_tox_env_observed = c(29, 27, 33, 5, 0)
 )
-plot_effect(model)
-plot_stress(model)
+plot_effect(model, show_legend = TRUE)
+plot_stress(model, show_legend = TRUE)
 
 # Plot all curves:
 plot_effect(model, which = "all")
 plot_stress(model, which = "all")
 
 # Plot only some selected curves:
-plot_effect(model, which = c("effect_tox", "effect_tox_env"))
-plot_stress(model, which = c("stress_tox", "stress_tox_env"))
+plot_effect(model, which = c("effect_tox_sys", "effect_tox_env_sys"))
+plot_stress(model, which = c("sys_tox", "sys_tox_env"))
+
+# Plot only the observed values:
+plot_effect(model, which = c("effect_tox_observed", "effect_tox_env_observed"))
+plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed"))
 }
diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R
index 01dcd29941a74c9616a370f3fd537601a2665b6f..d04e16f80ec3de4c862e7cde047edc5bfabda006 100644
--- a/tests/testthat/test-ecxsys.R
+++ b/tests/testthat/test-ecxsys.R
@@ -105,7 +105,7 @@ test_that("the discrete results have not changed", {
         tolerance = 1e-4
     )
     expect_equal(
-        mod$sys_tox_not_fitted,
+        mod$sys_tox_observed,
         c(0.2541156, 0.2324304, 0, 0, 0),
         tolerance = 1e-4
     )
@@ -131,7 +131,7 @@ test_that("the discrete results have not changed", {
         tolerance = 1e-4
     )
     expect_equal(
-        mod$sys_tox_env_not_fitted,
+        mod$sys_tox_env_observed,
         c(0.2566005, 0.1753250, 0, 0, 0),
         tolerance = 1e-4
     )