diff --git a/R/data.R b/R/data.R
index 9c5569646..48548ffed 100644
--- a/R/data.R
+++ b/R/data.R
@@ -217,5 +217,7 @@
 #' A named list with functions:
 #' - "wis" = [wis()]
 #' - "bias" = [bias_quantile()]
+#' - "coverage_50" = \(...) {run_safely(..., range = 50, fun = interval_coverage_quantile)} #nolint
+#' - "coverage_90" = \(...) {run_safely(..., range = 90, fun = interval_coverage_quantile)} #nolint
 #' @keywords info
-"metrics_sample"
+"metrics_quantile"
diff --git a/R/metrics-quantile.R b/R/metrics-quantile.R
index 5f45426ab..4f3774a96 100644
--- a/R/metrics-quantile.R
+++ b/R/metrics-quantile.R
@@ -149,6 +149,7 @@ interval_coverage_quantile <- function(observed, predicted, quantile, range = 50
 #' which predictions were made. If this does not contain the median (0.5) then
 #' the median is imputed as being the mean of the two innermost quantiles.
 #' @inheritParams bias_range
+#' @param na.rm logical. Should missing values be removed?
 #' @return scalar with the quantile bias for a single quantile prediction
 #' @author Nikos Bosse \email{nikosbosse@@gmail.com}
 #' @examples
@@ -167,69 +168,84 @@ interval_coverage_quantile <- function(observed, predicted, quantile, range = 50
 #' bias_quantile(observed, predicted, quantile)
 #' @export
 #' @keywords metric
-
-bias_quantile <- function(observed, predicted, quantile) {
-
+bias_quantile <- function(observed, predicted, quantile, na.rm = TRUE) {
   assert_input_quantile(observed, predicted, quantile)
-
   n <- length(observed)
   N <- length(quantile)
-
-  dt <- data.table(
-    forecast_id = rep(1:n, each = N),
-    observed = rep(observed, each = N),
-    predicted = as.vector(t(predicted)),
-    quantile = quantile
-  )[order(forecast_id, quantile)]
-
-  dt <- dt[, .(bias = bias_quantile_single(.SD)),   by = forecast_id]
-
-  return(dt$bias)
+  if (is.null(dim(predicted))) {
+    dim(predicted) <- c(n, N)
+  }
+  bias <- sapply(1:n, function(i) {
+    bias_quantile_single_vector(observed[i], predicted[i,], quantile, na.rm)
+  })
+  return(bias)
 }
 
 
-bias_quantile_single <- function(dt) {
+#' Compute Bias for a Single Vector of Quantile Predictions
+#' @description Internal function to compute bias for a single observed value,
+#' a vector of predicted values and a vector of quantiles.
+#' @param observed scalar with the observed value
+#' @param predicted vector of length N corresponding to the number of quantiles
+#' that holds predictions
+#' @param quantile vector of corresponding size N with the quantile levels for
+#' which predictions were made. If this does not contain the median (0.5) then
+#' the median is imputed as being the mean of the two innermost quantiles.
+#' @inheritParams bias_quantile
+#' @return scalar with the quantile bias for a single quantile prediction
+bias_quantile_single_vector <- function(observed, predicted, quantile, na.rm) {
 
-  dt <- dt[!is.na(quantile) & !is.na(predicted)]
+  assert_number(observed)
+  # other checks should have happend before
 
-  observed <- unique(dt$observed)
+  predicted_has_NAs <- anyNA(predicted)
+  quantile_has_NAs <- anyNA(quantile)
 
-  if (nrow(dt) == 0) {
-    return(NA_real_)
+  if(any(predicted_has_NAs, quantile_has_NAs)) {
+    if (!na.rm) {
+      return(NA_real_)
+    } else {
+      quantile <- quantile[!is.na(predicted)]
+      predicted <- predicted[!is.na(predicted)]
+      predicted <- predicted[!is.na(quantile)]
+      quantile <- quantile[!is.na(quantile)]
+    }
   }
 
-  if (!all(diff(dt$predicted) >= 0)) {
+  order <- order(quantile)
+  predicted <- predicted[order]
+  if (!all(diff(predicted) >= 0)) {
     stop("Predictions must not be decreasing with increasing quantile level")
   }
 
-  if (0.5 %in% dt$quantile) {
-    median_prediction <- dt[quantile == 0.5]$predicted
+  if (0.5 %in% quantile) {
+    median_prediction <- predicted[quantile == 0.5]
   } else {
-    # if median is not available, compute as mean of two innermost quantiles
+    # if median is not available, compute as mean of two innermost quantile
     message(
-      "Median not available, computing as mean of two innermost quantiles",
+      "Median not available, computing as mean of two innermost quantile",
       " in order to compute bias."
     )
     median_prediction <-
-      0.5 * dt[quantile == max(quantile[quantile < 0.5])]$predicted +
-      0.5 * dt[quantile == min(quantile[quantile > 0.5])]$predicted
+      0.5 * predicted[quantile == max(quantile[quantile < 0.5])] +
+      0.5 * predicted[quantile == min(quantile[quantile > 0.5])]
   }
 
   if (observed == median_prediction) {
     bias <- 0
     return(bias)
   } else if (observed < median_prediction) {
-    if (observed < min(dt$predicted)) {
+    if (observed < min(predicted)) {
       bias <- 1
     } else {
-      q <- max(dt[predicted <= observed]$quantile)
+      q <- max(quantile[predicted <= observed])
       bias <- 1 - 2 * q
     }
   } else if (observed > median_prediction) {
-    if (observed > max(dt$predicted)) {
+    if (observed > max(predicted)) {
       bias <- -1
     } else {
-      q <- min(dt[predicted >= observed]$quantile)
+      q <- min(quantile[predicted >= observed])
       bias <- 1 - 2 * q
     }
   }
diff --git a/data/metrics_quantile.rda b/data/metrics_quantile.rda
index c0217df2c..f1e231d77 100644
Binary files a/data/metrics_quantile.rda and b/data/metrics_quantile.rda differ
diff --git a/man/bias_quantile.Rd b/man/bias_quantile.Rd
index 4cbeb0338..a5f4b6492 100644
--- a/man/bias_quantile.Rd
+++ b/man/bias_quantile.Rd
@@ -4,7 +4,7 @@
 \alias{bias_quantile}
 \title{Determines Bias of Quantile Forecasts}
 \usage{
-bias_quantile(observed, predicted, quantile)
+bias_quantile(observed, predicted, quantile, na.rm = TRUE)
 }
 \arguments{
 \item{observed}{a single observed value}
@@ -15,6 +15,8 @@ that holds predictions}
 \item{quantile}{vector of corresponding size with the quantile levels for
 which predictions were made. If this does not contain the median (0.5) then
 the median is imputed as being the mean of the two innermost quantiles.}
+
+\item{na.rm}{logical. Should missing values be removed?}
 }
 \value{
 scalar with the quantile bias for a single quantile prediction
diff --git a/man/bias_quantile_single_vector.Rd b/man/bias_quantile_single_vector.Rd
new file mode 100644
index 000000000..26ac893b5
--- /dev/null
+++ b/man/bias_quantile_single_vector.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/metrics-quantile.R
+\name{bias_quantile_single_vector}
+\alias{bias_quantile_single_vector}
+\title{Compute Bias for a Single Vector of Quantile Predictions}
+\usage{
+bias_quantile_single_vector(observed, predicted, quantile, na.rm)
+}
+\arguments{
+\item{observed}{scalar with the observed value}
+
+\item{predicted}{vector of length N corresponding to the number of quantiles
+that holds predictions}
+
+\item{quantile}{vector of corresponding size N with the quantile levels for
+which predictions were made. If this does not contain the median (0.5) then
+the median is imputed as being the mean of the two innermost quantiles.}
+
+\item{na.rm}{logical. Should missing values be removed?}
+}
+\value{
+scalar with the quantile bias for a single quantile prediction
+}
+\description{
+Internal function to compute bias for a single observed value,
+a vector of predicted values and a vector of quantiles.
+}
diff --git a/man/metrics_quantile.Rd b/man/metrics_quantile.Rd
new file mode 100644
index 000000000..e96dcc836
--- /dev/null
+++ b/man/metrics_quantile.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{metrics_quantile}
+\alias{metrics_quantile}
+\title{Default metrics for quantile-based forecasts.}
+\format{
+An object of class \code{list} of length 4.
+}
+\usage{
+metrics_quantile
+}
+\description{
+A named list with functions:
+\itemize{
+\item "wis" = \code{\link[=wis]{wis()}}
+\item "bias" = \code{\link[=bias_quantile]{bias_quantile()}}
+\item "coverage_50" = \(...) {run_safely(..., range = 50, fun = interval_coverage_quantile)} #nolint
+\item "coverage_90" = \(...) {run_safely(..., range = 90, fun = interval_coverage_quantile)} #nolint
+}
+}
+\keyword{info}
diff --git a/man/metrics_sample.Rd b/man/metrics_sample.Rd
index 90aa265f7..5231f4ae7 100644
--- a/man/metrics_sample.Rd
+++ b/man/metrics_sample.Rd
@@ -5,13 +5,9 @@
 \alias{metrics_sample}
 \title{Default metrics for sample-based forecasts.}
 \format{
-An object of class \code{list} of length 7.
-
 An object of class \code{list} of length 7.
 }
 \usage{
-metrics_sample
-
 metrics_sample
 }
 \description{
@@ -26,11 +22,5 @@ A named list with functions:
 \item "ae_median" = \code{\link[=ae_median_sample]{ae_median_sample()}}
 \item "se_mean" = \code{\link[=se_mean_sample]{se_mean_sample()}}
 }
-
-A named list with functions:
-\itemize{
-\item "wis" = \code{\link[=wis]{wis()}}
-\item "bias" = \code{\link[=bias_quantile]{bias_quantile()}}
-}
 }
 \keyword{info}
diff --git a/tests/testthat/test-bias.R b/tests/testthat/test-bias.R
index b1e59fd6d..42eefa34f 100644
--- a/tests/testthat/test-bias.R
+++ b/tests/testthat/test-bias.R
@@ -129,6 +129,10 @@ test_that("bias_quantile() handles NA values", {
     bias_quantile(observed = 2, predicted, quantiles),
     -1
   )
+  expect_equal(
+    bias_quantile(observed = 2, predicted, quantiles, na.rm = FALSE),
+    NA_real_
+  )
 })
 
 test_that("bias_quantile() errors if no predictions", {