diff --git a/R/metrics-quantile.R b/R/metrics-quantile.R
index 3d9c76b9..5b57c382 100644
--- a/R/metrics-quantile.R
+++ b/R/metrics-quantile.R
@@ -275,16 +275,14 @@ interval_coverage <- function(observed, predicted,
   ) / 100
   if (!all(necessary_quantiles %in% quantile_level)) {
     #nolint start: keyword_quote_linter object_usage_linter
-    cli_warn(
+    cli_abort(
       c(
         "!" = "To compute the interval coverage for an interval range of
         {.val {interval_range}%}, the {.val {necessary_quantiles}} quantiles
-        are required.",
-        "i" = "Returning {.val {NA}}."
+        are required"
       )
     )
     #nolint end
-    return(NA)
   }
   r <- interval_range
   reformatted <- quantile_to_interval(observed, predicted, quantile_level)
@@ -332,7 +330,10 @@ interval_coverage <- function(observed, predicted,
 #' Bias can assume values between -1 and 1 and is 0 ideally (i.e. unbiased).
 #'
 #' Note that if the given quantiles do not contain the median, the median is
-#' imputed as the mean of the two innermost quantiles.
+#' imputed as a linear interpolation of the two innermost quantiles. If the
+#' median is not available and cannot be imputed, an error will be thrown.
+#' Note that in order to compute bias, quantiles must be non-decreasing with
+#' increasing quantile levels.
 #'
 #' For a large enough number of quantiles, the
 #' percentile rank will equal the proportion of predictive samples below the
@@ -357,6 +358,11 @@ interval_coverage <- function(observed, predicted,
 #' bias_quantile(observed, predicted, quantile_level)
 bias_quantile <- function(observed, predicted, quantile_level, na.rm = TRUE) {
   assert_input_quantile(observed, predicted, quantile_level)
+  # for bias quantile to work, at least one quantile level has to be <= 0.5
+  # and at least one >= 0.5
+  assert_vector(quantile_level[quantile_level <= 0.5], min.len = 1)
+  assert_vector(quantile_level[quantile_level >= 0.5], min.len = 1)
+
   n <- length(observed)
   N <- length(quantile_level)
   if (is.null(dim(predicted))) {
@@ -506,15 +512,13 @@ ae_median_quantile <- function(observed, predicted, quantile_level) {
   assert_input_quantile(observed, predicted, quantile_level)
   if (!any(quantile_level == 0.5)) {
     #nolint start: keyword_quote_linter
-    cli_warn(
+    cli_abort(
       c(
-        "x" = "In order to compute the absolute error of the median,
-        {.val 0.5} must be among the quantiles given.",
-        "i" = "Returning {.val NA}."
+        "!" = "In order to compute the absolute error of the median,
+        {.val 0.5} must be among the quantiles given"
       )
     )
     #nolint end
-    return(NA_real_)
   }
   if (is.null(dim(predicted))) {
     predicted <- matrix(predicted, nrow = 1)
diff --git a/man/bias_quantile.Rd b/man/bias_quantile.Rd
index fd45476d..428ba3af 100644
--- a/man/bias_quantile.Rd
+++ b/man/bias_quantile.Rd
@@ -64,7 +64,10 @@ out)
 Bias can assume values between -1 and 1 and is 0 ideally (i.e. unbiased).
 
 Note that if the given quantiles do not contain the median, the median is
-imputed as the mean of the two innermost quantiles.
+imputed as a linear interpolation of the two innermost quantiles. If the
+median is not available and cannot be imputed, an error will be thrown.
+Note that in order to compute bias, quantiles must be non-decreasing with
+increasing quantile levels.
 
 For a large enough number of quantiles, the
 percentile rank will equal the proportion of predictive samples below the
diff --git a/tests/testthat/test-metrics-quantile.R b/tests/testthat/test-metrics-quantile.R
index 16b9bfe1..0fbf8a5d 100644
--- a/tests/testthat/test-metrics-quantile.R
+++ b/tests/testthat/test-metrics-quantile.R
@@ -659,7 +659,7 @@ test_that("interval_coverage rejects wrong inputs", {
 test_that("interval_coverage_quantile throws a warning when a required quantile is not available", {
   dropped_quantile_pred <- predicted[, -4]
   dropped_quantiles <- quantile_level[-4]
-  expect_warning(
+  expect_error(
     interval_coverage(
       observed, dropped_quantile_pred, dropped_quantiles, interval_range = 50
     ),
@@ -858,6 +858,17 @@ test_that("bias_quantile() works with point forecasts", {
 })
 
 
+test_that("bias_quantile() handles cases where median is not available", {
+  predicted <- c(1, 10)
+  observed <- 15
+  quantile_level <- c(0.2, 0.4)
+
+  expect_error(
+    bias_quantile(observed, predicted, quantile_level),
+    "Assertion on 'quantile_level\\[quantile_leve\\l >= 0.5]' failed: Must have length >= 1, but has length 0."
+  )
+})
+
 # `interpolate_median` ======================================================= #
 test_that("interpolation in `interpolate_median` works", {
   predicted <- c(1, 10)
@@ -910,16 +921,10 @@ test_that("ae_median_quantile() works as_expected", {
   )
 
   # test that we get a warning if there are inputs without a 0.5 quantile
-  expect_warning(
-    expect_equal(
-      ae_median_quantile(observed, predicted_values, quantile_level = 0.6),
-      NA_real_
-    ),
-    'In order to compute the absolute error of the median, "0.5" must be among the quantiles given.'
+  expect_error(
+    ae_median_quantile(observed, predicted_values, quantile_level = 0.6),
+    'In order to compute the absolute error of the median, '
   )
-
-
-  dim(1:10)
 })
 
 
diff --git a/tests/testthat/test-score.R b/tests/testthat/test-score.R
index fb25f004..3a2f3594 100644
--- a/tests/testthat/test-score.R
+++ b/tests/testthat/test-score.R
@@ -304,6 +304,56 @@ test_that("score.forecast_quantile() works as expected in edge cases", {
 })
 
 
+
+test_that("score() works even if only some quantiles are missing", {
+
+  # only the median is there
+  onlymedian <- example_quantile[quantile_level == 0.5]
+  expect_no_condition(
+    score(onlymedian, metrics = get_metrics(
+      example_quantile,
+      exclude = c("interval_coverage_50", "interval_coverage_90")
+    ))
+  )
+
+
+  # asymmetric intervals
+  asymm <- example_quantile[!quantile_level > 0.6]
+  expect_warning(
+    expect_warning(
+      score_a <- score(asymm) %>% summarise_scores(by = "model"),
+      "Computation for `interval_coverage_50` failed."
+    ),
+    "Computation for `interval_coverage_90` failed."
+  )
+
+  # check that the result is equal to a case where we discard the entire
+  # interval in terms of WIS
+  inner <- example_quantile[quantile_level %in% c(0.4, 0.45, 0.5, 0.55, 0.6)]
+  score_b <- score(inner, get_metrics(
+    inner, exclude = c("interval_coverage_50", "interval_coverage_90")
+  )) %>%
+    summarise_scores(by = "model")
+  expect_equal(
+    score_a$wis,
+    score_b$wis
+  )
+
+  # median is not there, but only in a single model
+  test <- data.table::copy(example_quantile)
+  test_no_median <- test[model == "epiforecasts-EpiNow2" & !(quantile_level %in% c(0.5)), ]
+  test <- rbind(test[model != "epiforecasts-EpiNow2"], test_no_median)
+
+  test <- suppressWarnings(as_forecast_quantile(test))
+  expect_message(
+    expect_warning(
+      score(test),
+      "Computation for `ae_median` failed."
+    ),
+    "interpolating median from the two innermost quantiles"
+  )
+})
+
 # test integer and continuous case ---------------------------------------------
 test_that("function produces output for a continuous format case", {