diff --git a/NEWS.md b/NEWS.md
index 4426786c1efcf9bb0ddf17e435908e3c5fa9448d..0d9f1d7f2aa6f788d63bb88effdd2fd0a46d79d0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,7 @@
 * Renamed `predict_mixture()` to `multi_tox()`.
 * The argument `proportion_ca` in the mixture model `multi_tox()` was renamed and its value reversed. It is now called `sa_contribution` and specifies the proportion of stress addition in the calculation of toxicant stress. To convert your code from the old version use this equation: sa_contribution = 1 - proportion_ca.
 * Renamed `stress_tox_sam` to `stress_tox_sa` in the output of `multi_tox()`.
+* Fixed a bug where `plot_stress()` with argument `which = NULL` would result in an error. Now it correctly draws the axes without data.
 
 # stressaddition 2.7.0
 
diff --git a/R/plot_stress.R b/R/plot_stress.R
index 9e661d615425da91d66d9ca753e13fbeaea75398..8bbd4b2f8b1d08eb28fca29ab137e8969d645f60 100644
--- a/R/plot_stress.R
+++ b/R/plot_stress.R
@@ -58,8 +58,19 @@ plot_stress <- function(model,
         model$args$concentration[-1]
     )
 
-    curves_w <- curves[, which[!endsWith(which, "observed")]]
-    ymax <- if (NCOL(curves_w) == 0) 1 else max(curves_w, 1, na.rm = TRUE)
+    if (is.null(which)) {
+        ymax = 1
+    } else {
+        which_lines <- which[!endsWith(which, "observed")]
+        if (length(which_lines) == 0) {
+            ymax <- 1
+        } else {
+            ymax <- max(curves[, which_lines], 1, na.rm = TRUE)
+            # No need to include the observed stress in the call to max() no
+            # matter if those are in "which" or not because these vectors are
+            # clamped to [0, 1] anyway.
+        }
+    }
 
     plot(
         NA,