diff --git a/NEWS.md b/NEWS.md
index e1ee239eb..51104e477 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,6 @@
 # scoringutils (development version)
 
+- Fixed `bias_quantile()` returning wrong values when quantile levels were passed unsorted: predictions were reordered by quantile level but the quantile levels themselves were not, so predictions and levels became mispaired. Also fixed a crash ("argument is of length zero") when `na.rm = TRUE` removed all quantile levels on one side of the median; `bias_quantile()` now returns `NA` in this case, consistent with `na.rm = FALSE` (#1198).
 - Added `filter_scores()` and `impute_missing_scores()` for handling missing forecasts before summarisation. `filter_scores()` removes target combinations with insufficient model coverage, while `impute_missing_scores()` fills in missing scores using configurable strategies (worst, mean, NA, or reference model). Both use a strategy function pattern for extensibility. See `vignette("handling-missing-forecasts")` for details (#1122).
 - Added `plot_discrimination()` to visualise the discrimination ability of binary forecasts by plotting the distribution of predicted probabilities, stratified by the observed outcome. The function requires a `forecast_binary` object (created with `as_forecast_binary()`) (#942).
 - Fixed `summarise_scores()` producing a data.table with duplicate column names when the input `scores` object had no score columns (e.g. because every metric in `score()` warned and returned nothing). `summarise_scores()` now matches metric columns by exact name rather than regex partial match, and errors with a clear message when there is nothing to summarise (#1179).
diff --git a/R/metrics-quantile.R b/R/metrics-quantile.R
index e11cd00d9..03af7dc7f 100644
--- a/R/metrics-quantile.R
+++ b/R/metrics-quantile.R
@@ -475,10 +475,16 @@ bias_quantile_single_vector <- function(observed, predicted,
     predicted <- predicted[!is.na(predicted)]
     predicted <- predicted[!is.na(quantile_level)]
     quantile_level <- quantile_level[!is.na(quantile_level)]
+    # if NA removal leaves no quantile level on one side of the median, the
+    # median cannot be interpolated and no bias can be computed
+    if (!any(quantile_level <= 0.5) || !any(quantile_level >= 0.5)) {
+      return(NA_real_)
+    }
   }
 
   order <- order(quantile_level)
   predicted <- predicted[order]
+  quantile_level <- quantile_level[order]
   if (!all(diff(predicted) >= 0)) {
     cli_abort(
       c(
diff --git a/tests/testthat/test-metrics-quantile.R b/tests/testthat/test-metrics-quantile.R
index 354c48444..a60f445fe 100644
--- a/tests/testthat/test-metrics-quantile.R
+++ b/tests/testthat/test-metrics-quantile.R
@@ -751,6 +751,35 @@ test_that("bias_quantile() handles NA values", {
   )
 })
 
+test_that("bias_quantile() is invariant to the order of quantile levels", {
+  expect_equal( # nolint: expect_identical_linter
+    bias_quantile(observed = 1.5, c(3, 1, 2), c(0.75, 0.25, 0.5)),
+    0.5
+  )
+  expect_equal(
+    bias_quantile(observed = 1.5, c(3, 1, 2), c(0.75, 0.25, 0.5)),
+    bias_quantile(observed = 1.5, c(1, 2, 3), c(0.25, 0.5, 0.75))
+  )
+})
+
+test_that("bias_quantile() returns NA when na.rm removes one side of the median", {
+  expect_equal( # nolint: expect_identical_linter
+    bias_quantile(observed = 2, c(NA, NA, 3), c(0.25, 0.5, 0.75), na.rm = TRUE),
+    NA_real_
+  )
+  expect_equal( # nolint: expect_identical_linter
+    suppressMessages(
+      bias_quantile(observed = 2, c(NA, 3), c(0.25, 0.75), na.rm = TRUE)
+    ),
+    NA_real_
+  )
+  # if the median itself survives NA removal, the forecast can still be scored
+  expect_equal( # nolint: expect_identical_linter
+    bias_quantile(observed = 1, c(NA, 2, 3), c(0.25, 0.5, 0.75), na.rm = TRUE),
+    1
+  )
+})
+
 test_that("bias_quantile() errors if no predictions", {
   expect_error(
     bias_quantile(observed = 2, numeric(0), numeric(0)),