tests/testthat/test-metrics-quantile.R
In epiforecasts/scoringutils: Utilities for Scoring and Assessing Predictions
observed <- c(1, -15, 22)
predicted <- rbind(
c(-1, 0, 1, 2, 3),
c(-2, 1, 2, 2, 4),
c(-2, 0, 3, 3, 4)
)
quantile_level <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# covidHubUtils test:
y <- c(1, -15, 22)
forecast_quantiles_matrix <- rbind(
c(-1, 0, 1, 2, 3),
c(-2, 1, 2, 2, 4),
c(-2, 0, 3, 3, 4)
)
forecast_quantile_probs <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# ============================================================================ #
# Input handling ===============================================================
# ============================================================================ #
test_that("Input checking for quantile forecasts works", {
# everything correct
expect_no_condition(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level)
)
# quantile_level > 1
expect_error(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level + 1),
"Assertion on 'quantile_level' failed: Element 1 is not <= 1."
)
# quantile_level < 0
expect_error(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level - 1),
"Assertion on 'quantile_level' failed: Element 1 is not >= 0."
)
# 10 observations, but only 3 forecasts
expect_error(
scoringutils:::assert_input_quantile(1:10, predicted, quantile_level),
"Assertion on 'predicted' failed: Must have exactly 10 rows, but has 3 rows."
)
# observed value is a factor
expect_error(
scoringutils:::assert_input_quantile(factor(1:10), predicted, quantile_level),
"Assertion on 'observed' failed: Must be of type 'numeric', not 'factor'."
)
# observed is a single number and does not have the same length as predicted
# There seems to be an issue with the error message: there is one \n to many
# such that the test fails when executed alone, but works when executed
# together with others.
expect_error(
scoringutils:::assert_input_quantile(1, predicted, quantile_level),
"Assertion failed. One of the following must apply:\n * check_numeric_vector(predicted): Must be of type 'atomic vector',\n * not 'matrix'\n * check_matrix(predicted): Must have exactly 1 rows, but has 3 rows",
fixed = TRUE
)
# predicted is a vector
expect_error(
scoringutils:::assert_input_quantile(observed, as.vector(predicted), quantile_level),
"Assertion on 'predicted' failed: Must be of type 'matrix', not 'double'."
)
})
# ============================================================================ #
# wis ==========================================================================
# ============================================================================ #a
test_that("wis works standalone, median only", {
y <- c(1, -15, 22)
lower <- upper <- predicted_quantile <- c(1, 2, 3)
quantile_probs <- 0.5
actual <- interval_score(y,
lower = lower, upper = upper,
weigh = TRUE,
interval_range = 0
)
actual_wis <- wis(
observed = y,
predicted = matrix(predicted_quantile),
quantile_level = quantile_probs,
)
expected <- abs(y - lower)
expect_identical(actual, expected)
})
test_that("`wis()` works within score for median forecast", {
test_data <- data.frame(
observed = c(1, -15, 22),
predicted = 1:3,
quantile_level = rep(c(0.5), each = 3),
model = "model1",
date = 1:3
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("ae_median", "wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
expect_equal(eval$ae_median, eval$wis)
})
test_that("`wis()` equals `interval_score()`, 1 interval only", {
y <- c(1, -15, 22)
lower <- c(0, 1, 0)
upper <- c(2, 2, 3)
quantile_probs <- c(0.25, 0.75)
predicted <- matrix(c(lower, upper), ncol = 2)
alpha <- 0.5
expected <- (upper - lower) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
actual <- interval_score(
y,
lower = lower, upper = upper,
weigh = TRUE,
interval_range = 50
)
actual_wis <- wis(
observed = y,
predicted = predicted,
quantile_level = quantile_probs,
)
expect_identical(actual, expected)
expect_identical(actual_wis, expected)
})
test_that("wis() works within score for one interval", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 2),
quantile_level = rep(c(0.25, 0.75), each = 3),
predicted = c(c(0, 1, 0), c(2, 2, 3)),
model = c("model1"),
date = rep(1:3, times = 2)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
lower <- c(0, 1, 0)
upper <- c(2, 2, 3)
alpha <- 0.5
expected <- (upper - lower) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
expect_equal(expected, eval$wis)
})
test_that("`wis()` works 1 interval and median", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 3),
quantile_level = rep(c(0.25, 0.5, 0.75), each = 3),
predicted = c(c(0, 1, 0), c(1, 2, 3), c(2, 2, 3)),
model = c("model1"),
date = rep(1:3, times = 3)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
y <- c(1, -15, 22)
quantile <- rbind(c(0, 1, 2), c(1, 2, 2), c(0, 3, 3))
quantile_probs <- c(0.25, 0.5, 0.75)
alpha <- 0.5
expected <- 0.5 * (
abs(y - quantile[, 2]) +
(quantile[, 3] - quantile[, 1]) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = quantile,
quantile_level = quantile_probs,
count_median_twice = TRUE
)
expect_identical(eval$wis, expected)
expect_identical(actual_wis, expected)
})
test_that("wis works, 2 intervals and median", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 5),
quantile_level = rep(c(0.1, 0.25, 0.5, 0.75, 0.9), each = 3),
predicted = c(
c(-1, -2, -2), c(0, 1, 0), c(1, 2, 3),
c(2, 2, 3), c(3, 4, 4)
),
model = c("model1"),
date = rep(1:3, times = 5)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
quantile <- forecast_quantiles_matrix
quantile_probs <- c(0.1, 0.25, 0.5, 0.75, 0.9)
alpha1 <- 0.2
alpha2 <- 0.5
expected <- (1 / 3) * (
abs(y - quantile[, 3]) +
(quantile[, 5] - quantile[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(quantile[, 4] - quantile[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = quantile,
quantile_level = quantile_probs,
count_median_twice = TRUE
)
expect_equal(
as.numeric(eval$wis),
as.numeric(expected)
)
expect_identical(actual_wis, expected)
})
# additional tests from the covidhubutils repo
test_that("wis is correct, median only - test corresponds to covidHubUtils", {
forecast_quantiles_matrix <- forecast_quantiles_matrix[, 3, drop = FALSE]
forecast_quantile_probs <- forecast_quantile_probs[3]
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted)
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(
as_forecast(data_formatted), metrics = metrics
)
expected <- abs(y - forecast_quantiles_matrix[, 1])
actual_wis <- wis(
observed = y,
predicted = matrix(forecast_quantiles_matrix),
quantile_level = 0.5,
count_median_twice = FALSE
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
})
test_that("wis is correct, 1 interval only - test corresponds to covidHubUtils", {
forecast_quantiles_matrix <- forecast_quantiles_matrix[, c(1, 5), drop = FALSE]
forecast_quantile_probs <- forecast_quantile_probs[c(1, 5)]
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
eval <- suppressMessages(score(data_formatted, metrics = metrics))
eval <- summarise_scores(eval,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
alpha1 <- 0.2
expected <- (forecast_quantiles_matrix[, 2] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) +
c(0, (-2) - (-15), 22 - 4)
actual_wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.9)
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
})
test_that("wis is correct, 2 intervals and median - test corresponds to covidHubUtils", {
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted) %>%
as_forecast()
eval <- score(data_formatted, metrics = metrics_no_cov)
eval <- summarise_scores(eval,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
alpha1 <- 0.2
alpha2 <- 0.5
expected <- (1 / 2.5) * (
0.5 * abs(y - forecast_quantiles_matrix[, 3]) +
(forecast_quantiles_matrix[, 5] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(forecast_quantiles_matrix[, 4] - forecast_quantiles_matrix[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.25, 0.5, 0.75, 0.9),
count_median_twice = FALSE
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
# check whether `count_median_twice = TRUE` also works
expected2 <- (1 / 3) * (
abs(y - forecast_quantiles_matrix[, 3]) +
(forecast_quantiles_matrix[, 5] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(forecast_quantiles_matrix[, 4] - forecast_quantiles_matrix[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis2 <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.25, 0.5, 0.75, 0.9),
count_median_twice = TRUE
)
metrics <- metrics_quantile() %>%
select_metrics("wis")
metrics$wis <- customise_metric(wis, count_median_twice = TRUE)
eval2 <- eval <- score(data_formatted, metrics = metrics)
eval2 <- summarise_scores(eval2,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
expect_equal(eval2$wis, expected2)
expect_equal(actual_wis2, expected2)
})
test_that("Quantlie score and interval score yield the same result, weigh = FALSE", {
# calling quantile_score and wis should return the same result if all
# quantiles form central prediction intervals
expect_equal(
quantile_score(observed, predicted, quantile_level),
wis(observed, predicted, quantile_level)
)
observed <- rnorm(10, mean = 1:10)
alphas <- c(0.1, 0.5, 0.9)
for (alpha in alphas) {
lower <- qnorm(alpha / 2, rnorm(10, mean = 1:10))
upper <- qnorm((1 - alpha / 2), rnorm(10, mean = 11:20))
w <- FALSE
is <- interval_score(
observed = observed,
lower = lower,
upper = upper,
interval_range = (1 - alpha) * 100,
weigh = w
)
wis <- wis(
observed = observed,
predicted = cbind(lower, upper),
quantile_level = c(alpha / 2, 1 - alpha / 2),
count_median_twice = FALSE,
weigh = w
)
qs_lower <- quantile_score(observed,
predicted = matrix(lower),
quantile_level = alpha / 2,
weigh = w
)
qs_upper <- quantile_score(observed,
predicted = matrix(upper),
quantile_level = 1 - alpha / 2,
weigh = w
)
expect_equal((qs_lower + qs_upper) / 2, is)
expect_equal(wis, is)
}
})
# ============================================================================ #
# Quantile score ============================================================= #
# ============================================================================ #
test_that("Quantlie score and interval score yield the same result, weigh = TRUE", {
observed <- rnorm(10, mean = 1:10)
alphas <- c(0.1, 0.5, 0.9)
for (alpha in alphas) {
lower <- qnorm(alpha / 2, rnorm(10, mean = 1:10))
upper <- qnorm((1 - alpha / 2), rnorm(10, mean = 11:20))
w <- TRUE
is <- interval_score(
observed = observed,
lower = lower,
upper = upper,
interval_range = (1 - alpha) * 100,
weigh = w
)
wis <- wis(
observed = observed,
predicted = cbind(lower, upper),
quantile_level = c(alpha / 2, 1 - alpha / 2),
count_median_twice = FALSE,
weigh = w
)
qs_lower <- quantile_score(observed,
predicted = matrix(lower),
quantile_level = alpha / 2,
weigh = w
)
qs_upper <- quantile_score(observed,
predicted = matrix(upper),
quantile_level = 1 - alpha / 2,
weigh = w
)
expect_equal((qs_lower + qs_upper) / 2, is)
expect_equal(wis, is)
}
})
test_that("wis works with separate results", {
wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = forecast_quantile_probs,
separate_results = TRUE
)
expect_equal(wis$wis, wis$dispersion + wis$overprediction + wis$underprediction)
})
# ============================================================================ #
# overprediction, underprediction, dispersion ================================ #
# ============================================================================ #
test_that("wis is the sum of overprediction, underprediction, dispersion", {
wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = forecast_quantile_probs
)
d <- dispersion(y, forecast_quantiles_matrix, forecast_quantile_probs)
o <- overprediction(y, forecast_quantiles_matrix, forecast_quantile_probs)
u <- underprediction(y, forecast_quantiles_matrix, forecast_quantile_probs)
expect_equal(wis, d + o + u)
})
# ============================================================================ #
# `interval_coverage` =============================================== #
# ============================================================================ #
test_that("interval_coverage works", {
expect_equal(
interval_coverage(observed, predicted, quantile_level, interval_range = 50),
c(TRUE, FALSE, FALSE)
)
})
test_that("interval_coverage rejects wrong inputs", {
expect_error(
interval_coverage(observed, predicted, quantile_level, interval_range = c(50, 0)),
"Assertion on 'interval_range' failed: Must have length 1."
)
})
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(
interval_coverage(
observed, dropped_quantile_pred, dropped_quantiles, interval_range = 50
),
paste(
"To compute the interval coverage for an interval range of \"50%\",",
"the 0.25 and 0.75 quantiles are required"
)
)
})
# ============================================================================ #
# `interval_coverage_deviation` ============================================== #
# ============================================================================ #
test_that("interval_coverage_deviation works", {
existing_ranges <- unique(get_range_from_quantile(quantile_level))
expect_equal(existing_ranges, c(80, 50, 0))
cov_50 <- interval_coverage(observed, predicted, quantile_level, interval_range = c(50))
cov_80 <- interval_coverage(observed, predicted, quantile_level, interval_range = c(80))
manual <- 0.5 * (cov_50 - 0.5) + 0.5 * (cov_80 - 0.8)
expect_equal(
interval_coverage_deviation(observed, predicted, quantile_level),
manual
)
expect_warning(
interval_coverage_deviation(
observed, predicted, c(quantile_level[-4], 0.76)
),
"all quantiles must form central symmetric prediction intervals"
)
})
# ============================================================================ #
# `bias_quantile` ============================================================ #
# ============================================================================ #
test_that("bias_quantile() works as expected", {
predicted <- c(1, 2, 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
0
)
predicted <- c(0, 1, 2)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
-0.8
)
predicted <- c(
705.500, 1127.000, 4006.250, 4341.500, 4709.000, 4821.996,
5340.500, 5451.000, 5703.500, 6087.014, 6329.500, 6341.000,
6352.500, 6594.986, 6978.500, 7231.000, 7341.500, 7860.004,
7973.000, 8340.500, 8675.750, 11555.000, 11976.500
)
quantile_level <- c(0.01, 0.025, seq(0.05, 0.95, 0.05), 0.975, 0.99)
observed <- 8062
expect_equal(bias_quantile(observed, predicted, quantile_level), -0.8)
})
test_that("bias_quantile handles matrix input", {
observed <- seq(10, 0, length.out = 4)
predicted <- matrix(1:12, ncol = 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed, predicted, quantiles),
c(-1.0, -0.8, 0.8, 1.0)
)
})
test_that("bias_quantile() handles vector that is too long", {
predicted <- c(NA, 1, 2, 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_error(
bias_quantile(observed = 2, predicted, quantiles),
"Assertion on 'quantile_level' failed: Must have length 4, but has length 3."
)
})
test_that("bias_quantile() handles NA values", {
predicted <- c(NA, 1, 2)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
-0.8
)
predicted <- c(0, 1, 2)
quantiles <- c(0.1, 0.5, NA)
expect_equal(
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", {
expect_error(
bias_quantile(observed = 2, numeric(0), numeric(0)),
"Assertion on 'quantile_level' failed: Must have length >= 1, but has length 0"
)
})
test_that("bias_quantile() returns correct bias if value below the median", {
predicted <- c(1, 2, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 1, predicted, quantiles), 0.8)
)
})
test_that("bias_quantile() returns correct bias if value above median", {
predicted <- c(1, 2, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 5, predicted, quantiles), -0.8)
)
})
test_that("bias_quantile() returns correct bias if value at the median", {
predicted <- c(1, 2, 3, 4)
quantiles <- c(0.1, 0.3, 0.5, 0.7)
expect_equal(bias_quantile(observed = 3, predicted, quantiles), 0)
})
test_that("bias_quantile() returns 1 if true value below min prediction", {
predicted <- c(2, 3, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 1, predicted, quantiles), 1)
)
})
test_that("bias_quantile() returns -1 if true value above max prediction", {
predicted <- c(1, 2, 3, 4)
quantiles <- c(0.1, 0.3, 0.5, 0.7)
expect_equal(bias_quantile(observed = 6, predicted, quantiles), -1)
})
test_that("bias_quantile(): quantiles must be between 0 and 1", {
predicted <- 1:4
# Failing example
quantiles <- c(-0.1, 0.3, 0.5, 0.8)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Assertion on 'quantile_level' failed: Element 1 is not >= 0."
)
# Passing counter example
quantiles <- c(0.1, 0.3, 0.5, 0.8)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile(): quantiles must be increasing", {
predicted <- 1:4
# Failing example
quantiles <- c(0.8, 0.3, 0.5, 0.9)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Predictions must not be decreasing with increasing quantile level"
)
# Passing counter example
quantiles <- c(0.3, 0.5, 0.8, 0.9)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile(): predictions must be increasing", {
predicted <- c(1, 2, 4, 3)
quantiles <- c(0.1, 0.3, 0.5, 0.9)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Predictions must not be decreasing with increasing quantile level"
)
expect_silent(bias_quantile( observed = 3, 1:4, quantiles))
})
test_that("bias_quantile(): quantile levels must be unique", {
predicted <- 1:4
# Failing example
quantiles <- c(0.3, 0.3, 0.5, 0.8)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Assertion on 'quantile_level' failed: Contains duplicated values, position 2."
)
# Passing example
quantiles <- c(0.3, 0.5, 0.8, 0.9)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile only produces one message", {
expect_message(
bias_quantile(observed, predicted[, -3], quantile_level[-3]),
"Median not available, computing bias as mean of the two innermost quantiles in order to compute bias."
)
})
test_that("bias_quantile() works with point forecasts", {
predicted <- 1
observed <- 1
quantile_level <- 0.5
expect_equal(bias_quantile(observed, predicted, quantile_level), 0)
})
# `interpolate_median` ======================================================= #
test_that("interpolation in `interpolate_median` works", {
predicted <- c(1, 10)
observed <- 15
quantile_level <- c(0.4, 0.6)
# median is missing, symmetric quantile levels given
expect_equal(interpolate_median(predicted, quantile_level), mean(predicted))
quantile_level <- c(0.1, 0.9)
expect_equal(interpolate_median(predicted, quantile_level), mean(predicted))
1 + 0.4/0.8 * (10 - 1)
# asymmetric quantile levels given
quantile_level <- c(0.3, 0.6)
expect_equal(interpolate_median(predicted, quantile_level), 1 + 0.2/0.3 * (10 - 1))
quantile_level <- c(0.2, 0.6)
expect_equal(interpolate_median(predicted, quantile_level), 1 + 0.3/0.4 * (10 - 1))
})
# ============================================================================ #
# `ae_median_quantile` ======================================================= #
# ============================================================================ #
test_that("ae_median_quantile() works as_expected", {
observed <- 1:30
predicted_values <- matrix(2:31)
expect_equal(
ae_median_quantile(observed, predicted_values, quantile_level = 0.5),
as.numeric(predicted_values - observed)
)
# using a vector as input for predicted does not work if there are several
# observed values
expect_error(
ae_median_quantile(observed, as.numeric(predicted_values), quantile_level = 0.5),
"Assertion on 'predicted' failed: Must be of type 'matrix', not 'double'."
)
# it does work if there is only a single observed value - in this case
# the vector gets treated as one forecast - which means that we need to
# supply multiple quantile levels now, because it assumes that single forecast
# is represented by several quantiles
expect_equal(
ae_median_quantile(observed[1], as.numeric(predicted_values)[1:3], quantile_level = c(0.1, 0.5, 0.7)),
2
)
# 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.'
)
dim(1:10)
})
epiforecasts/scoringutils documentation built on April 23, 2024, 4:56 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
observed <- c(1, -15, 22)
predicted <- rbind(
c(-1, 0, 1, 2, 3),
c(-2, 1, 2, 2, 4),
c(-2, 0, 3, 3, 4)
)
quantile_level <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# covidHubUtils test:
y <- c(1, -15, 22)
forecast_quantiles_matrix <- rbind(
c(-1, 0, 1, 2, 3),
c(-2, 1, 2, 2, 4),
c(-2, 0, 3, 3, 4)
)
forecast_quantile_probs <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# ============================================================================ #
# Input handling ===============================================================
# ============================================================================ #
test_that("Input checking for quantile forecasts works", {
# everything correct
expect_no_condition(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level)
)
# quantile_level > 1
expect_error(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level + 1),
"Assertion on 'quantile_level' failed: Element 1 is not <= 1."
)
# quantile_level < 0
expect_error(
scoringutils:::assert_input_quantile(observed, predicted, quantile_level - 1),
"Assertion on 'quantile_level' failed: Element 1 is not >= 0."
)
# 10 observations, but only 3 forecasts
expect_error(
scoringutils:::assert_input_quantile(1:10, predicted, quantile_level),
"Assertion on 'predicted' failed: Must have exactly 10 rows, but has 3 rows."
)
# observed value is a factor
expect_error(
scoringutils:::assert_input_quantile(factor(1:10), predicted, quantile_level),
"Assertion on 'observed' failed: Must be of type 'numeric', not 'factor'."
)
# observed is a single number and does not have the same length as predicted
# There seems to be an issue with the error message: there is one \n to many
# such that the test fails when executed alone, but works when executed
# together with others.
expect_error(
scoringutils:::assert_input_quantile(1, predicted, quantile_level),
"Assertion failed. One of the following must apply:\n * check_numeric_vector(predicted): Must be of type 'atomic vector',\n * not 'matrix'\n * check_matrix(predicted): Must have exactly 1 rows, but has 3 rows",
fixed = TRUE
)
# predicted is a vector
expect_error(
scoringutils:::assert_input_quantile(observed, as.vector(predicted), quantile_level),
"Assertion on 'predicted' failed: Must be of type 'matrix', not 'double'."
)
})
# ============================================================================ #
# wis ==========================================================================
# ============================================================================ #a
test_that("wis works standalone, median only", {
y <- c(1, -15, 22)
lower <- upper <- predicted_quantile <- c(1, 2, 3)
quantile_probs <- 0.5
actual <- interval_score(y,
lower = lower, upper = upper,
weigh = TRUE,
interval_range = 0
)
actual_wis <- wis(
observed = y,
predicted = matrix(predicted_quantile),
quantile_level = quantile_probs,
)
expected <- abs(y - lower)
expect_identical(actual, expected)
})
test_that("`wis()` works within score for median forecast", {
test_data <- data.frame(
observed = c(1, -15, 22),
predicted = 1:3,
quantile_level = rep(c(0.5), each = 3),
model = "model1",
date = 1:3
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("ae_median", "wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
expect_equal(eval$ae_median, eval$wis)
})
test_that("`wis()` equals `interval_score()`, 1 interval only", {
y <- c(1, -15, 22)
lower <- c(0, 1, 0)
upper <- c(2, 2, 3)
quantile_probs <- c(0.25, 0.75)
predicted <- matrix(c(lower, upper), ncol = 2)
alpha <- 0.5
expected <- (upper - lower) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
actual <- interval_score(
y,
lower = lower, upper = upper,
weigh = TRUE,
interval_range = 50
)
actual_wis <- wis(
observed = y,
predicted = predicted,
quantile_level = quantile_probs,
)
expect_identical(actual, expected)
expect_identical(actual_wis, expected)
})
test_that("wis() works within score for one interval", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 2),
quantile_level = rep(c(0.25, 0.75), each = 3),
predicted = c(c(0, 1, 0), c(2, 2, 3)),
model = c("model1"),
date = rep(1:3, times = 2)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
lower <- c(0, 1, 0)
upper <- c(2, 2, 3)
alpha <- 0.5
expected <- (upper - lower) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
expect_equal(expected, eval$wis)
})
test_that("`wis()` works 1 interval and median", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 3),
quantile_level = rep(c(0.25, 0.5, 0.75), each = 3),
predicted = c(c(0, 1, 0), c(1, 2, 3), c(2, 2, 3)),
model = c("model1"),
date = rep(1:3, times = 3)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
y <- c(1, -15, 22)
quantile <- rbind(c(0, 1, 2), c(1, 2, 2), c(0, 3, 3))
quantile_probs <- c(0.25, 0.5, 0.75)
alpha <- 0.5
expected <- 0.5 * (
abs(y - quantile[, 2]) +
(quantile[, 3] - quantile[, 1]) * (alpha / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = quantile,
quantile_level = quantile_probs,
count_median_twice = TRUE
)
expect_identical(eval$wis, expected)
expect_identical(actual_wis, expected)
})
test_that("wis works, 2 intervals and median", {
test_data <- data.frame(
observed = rep(c(1, -15, 22), times = 5),
quantile_level = rep(c(0.1, 0.25, 0.5, 0.75, 0.9), each = 3),
predicted = c(
c(-1, -2, -2), c(0, 1, 0), c(1, 2, 3),
c(2, 2, 3), c(3, 4, 4)
),
model = c("model1"),
date = rep(1:3, times = 5)
) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(test_data, metrics = metrics)
eval <- summarise_scores(eval, by = c("model", "date"))
quantile <- forecast_quantiles_matrix
quantile_probs <- c(0.1, 0.25, 0.5, 0.75, 0.9)
alpha1 <- 0.2
alpha2 <- 0.5
expected <- (1 / 3) * (
abs(y - quantile[, 3]) +
(quantile[, 5] - quantile[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(quantile[, 4] - quantile[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = quantile,
quantile_level = quantile_probs,
count_median_twice = TRUE
)
expect_equal(
as.numeric(eval$wis),
as.numeric(expected)
)
expect_identical(actual_wis, expected)
})
# additional tests from the covidhubutils repo
test_that("wis is correct, median only - test corresponds to covidHubUtils", {
forecast_quantiles_matrix <- forecast_quantiles_matrix[, 3, drop = FALSE]
forecast_quantile_probs <- forecast_quantile_probs[3]
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted)
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
metrics$wis <- customise_metric(metrics$wis, count_median_twice = TRUE)
eval <- score(
as_forecast(data_formatted), metrics = metrics
)
expected <- abs(y - forecast_quantiles_matrix[, 1])
actual_wis <- wis(
observed = y,
predicted = matrix(forecast_quantiles_matrix),
quantile_level = 0.5,
count_median_twice = FALSE
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
})
test_that("wis is correct, 1 interval only - test corresponds to covidHubUtils", {
forecast_quantiles_matrix <- forecast_quantiles_matrix[, c(1, 5), drop = FALSE]
forecast_quantile_probs <- forecast_quantile_probs[c(1, 5)]
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted) %>%
as_forecast()
metrics <- metrics_quantile() %>%
select_metrics(select = c("wis"))
eval <- suppressMessages(score(data_formatted, metrics = metrics))
eval <- summarise_scores(eval,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
alpha1 <- 0.2
expected <- (forecast_quantiles_matrix[, 2] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) +
c(0, (-2) - (-15), 22 - 4)
actual_wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.9)
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
})
test_that("wis is correct, 2 intervals and median - test corresponds to covidHubUtils", {
target_end_dates <- as.Date("2020-01-01") + c(7, 14, 7)
horizons <- c("1", "2", "1")
locations <- c("01", "01", "02")
target_variables <- rep("inc death", length(y))
forecast_target_end_dates <-
rep(target_end_dates, times = ncol(forecast_quantiles_matrix))
forecast_horizons <- rep(horizons, times = ncol(forecast_quantiles_matrix))
forecast_locations <- rep(locations, times = ncol(forecast_quantiles_matrix))
forecast_target_variables <-
rep(target_variables, times = ncol(forecast_quantiles_matrix))
forecast_quantile_probs <- rep(forecast_quantile_probs, each = length(y))
forecast_quantiles <- forecast_quantiles_matrix
dim(forecast_quantiles) <- prod(dim(forecast_quantiles))
test_truth <- data.frame(
model = rep("truth_source", length(y)),
target_variable = target_variables,
target_end_date = target_end_dates,
location = locations,
value = y,
stringsAsFactors = FALSE
)
n_forecasts <- length(forecast_quantiles)
test_forecasts <- data.frame(
model = rep("m1", n_forecasts),
forecast_date = rep(as.Date("2020-01-01"), n_forecasts),
location = forecast_locations,
horizon = forecast_horizons,
temporal_resolution = rep("wk", n_forecasts),
target_variable = forecast_target_variables,
target_end_date = forecast_target_end_dates,
type = rep("quantile", n_forecasts),
quantile_level = forecast_quantile_probs,
value = forecast_quantiles,
stringsAsFactors = FALSE
)
# make a version that conforms to scoringutils format
truth_formatted <- data.table::as.data.table(test_truth)
truth_formatted[, `:=`(model = NULL)]
data.table::setnames(truth_formatted, old = "value", new = "observed")
forecasts_formated <- data.table::as.data.table(test_forecasts)
data.table::setnames(forecasts_formated, old = "value", new = "predicted")
data_formatted <- merge(forecasts_formated, truth_formatted) %>%
as_forecast()
eval <- score(data_formatted, metrics = metrics_no_cov)
eval <- summarise_scores(eval,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
alpha1 <- 0.2
alpha2 <- 0.5
expected <- (1 / 2.5) * (
0.5 * abs(y - forecast_quantiles_matrix[, 3]) +
(forecast_quantiles_matrix[, 5] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(forecast_quantiles_matrix[, 4] - forecast_quantiles_matrix[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.25, 0.5, 0.75, 0.9),
count_median_twice = FALSE
)
expect_equal(eval$wis, expected)
expect_equal(actual_wis, expected)
# check whether `count_median_twice = TRUE` also works
expected2 <- (1 / 3) * (
abs(y - forecast_quantiles_matrix[, 3]) +
(forecast_quantiles_matrix[, 5] - forecast_quantiles_matrix[, 1]) * (alpha1 / 2) + c(0, (-2) - (-15), 22 - 4) +
(forecast_quantiles_matrix[, 4] - forecast_quantiles_matrix[, 2]) * (alpha2 / 2) + c(0, 1 - (-15), 22 - 3)
)
actual_wis2 <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = c(0.1, 0.25, 0.5, 0.75, 0.9),
count_median_twice = TRUE
)
metrics <- metrics_quantile() %>%
select_metrics("wis")
metrics$wis <- customise_metric(wis, count_median_twice = TRUE)
eval2 <- eval <- score(data_formatted, metrics = metrics)
eval2 <- summarise_scores(eval2,
by = c(
"model", "location", "target_variable",
"target_end_date", "forecast_date", "horizon"
)
)
expect_equal(eval2$wis, expected2)
expect_equal(actual_wis2, expected2)
})
test_that("Quantlie score and interval score yield the same result, weigh = FALSE", {
# calling quantile_score and wis should return the same result if all
# quantiles form central prediction intervals
expect_equal(
quantile_score(observed, predicted, quantile_level),
wis(observed, predicted, quantile_level)
)
observed <- rnorm(10, mean = 1:10)
alphas <- c(0.1, 0.5, 0.9)
for (alpha in alphas) {
lower <- qnorm(alpha / 2, rnorm(10, mean = 1:10))
upper <- qnorm((1 - alpha / 2), rnorm(10, mean = 11:20))
w <- FALSE
is <- interval_score(
observed = observed,
lower = lower,
upper = upper,
interval_range = (1 - alpha) * 100,
weigh = w
)
wis <- wis(
observed = observed,
predicted = cbind(lower, upper),
quantile_level = c(alpha / 2, 1 - alpha / 2),
count_median_twice = FALSE,
weigh = w
)
qs_lower <- quantile_score(observed,
predicted = matrix(lower),
quantile_level = alpha / 2,
weigh = w
)
qs_upper <- quantile_score(observed,
predicted = matrix(upper),
quantile_level = 1 - alpha / 2,
weigh = w
)
expect_equal((qs_lower + qs_upper) / 2, is)
expect_equal(wis, is)
}
})
# ============================================================================ #
# Quantile score ============================================================= #
# ============================================================================ #
test_that("Quantlie score and interval score yield the same result, weigh = TRUE", {
observed <- rnorm(10, mean = 1:10)
alphas <- c(0.1, 0.5, 0.9)
for (alpha in alphas) {
lower <- qnorm(alpha / 2, rnorm(10, mean = 1:10))
upper <- qnorm((1 - alpha / 2), rnorm(10, mean = 11:20))
w <- TRUE
is <- interval_score(
observed = observed,
lower = lower,
upper = upper,
interval_range = (1 - alpha) * 100,
weigh = w
)
wis <- wis(
observed = observed,
predicted = cbind(lower, upper),
quantile_level = c(alpha / 2, 1 - alpha / 2),
count_median_twice = FALSE,
weigh = w
)
qs_lower <- quantile_score(observed,
predicted = matrix(lower),
quantile_level = alpha / 2,
weigh = w
)
qs_upper <- quantile_score(observed,
predicted = matrix(upper),
quantile_level = 1 - alpha / 2,
weigh = w
)
expect_equal((qs_lower + qs_upper) / 2, is)
expect_equal(wis, is)
}
})
test_that("wis works with separate results", {
wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = forecast_quantile_probs,
separate_results = TRUE
)
expect_equal(wis$wis, wis$dispersion + wis$overprediction + wis$underprediction)
})
# ============================================================================ #
# overprediction, underprediction, dispersion ================================ #
# ============================================================================ #
test_that("wis is the sum of overprediction, underprediction, dispersion", {
wis <- wis(
observed = y,
predicted = forecast_quantiles_matrix,
quantile_level = forecast_quantile_probs
)
d <- dispersion(y, forecast_quantiles_matrix, forecast_quantile_probs)
o <- overprediction(y, forecast_quantiles_matrix, forecast_quantile_probs)
u <- underprediction(y, forecast_quantiles_matrix, forecast_quantile_probs)
expect_equal(wis, d + o + u)
})
# ============================================================================ #
# `interval_coverage` =============================================== #
# ============================================================================ #
test_that("interval_coverage works", {
expect_equal(
interval_coverage(observed, predicted, quantile_level, interval_range = 50),
c(TRUE, FALSE, FALSE)
)
})
test_that("interval_coverage rejects wrong inputs", {
expect_error(
interval_coverage(observed, predicted, quantile_level, interval_range = c(50, 0)),
"Assertion on 'interval_range' failed: Must have length 1."
)
})
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(
interval_coverage(
observed, dropped_quantile_pred, dropped_quantiles, interval_range = 50
),
paste(
"To compute the interval coverage for an interval range of \"50%\",",
"the 0.25 and 0.75 quantiles are required"
)
)
})
# ============================================================================ #
# `interval_coverage_deviation` ============================================== #
# ============================================================================ #
test_that("interval_coverage_deviation works", {
existing_ranges <- unique(get_range_from_quantile(quantile_level))
expect_equal(existing_ranges, c(80, 50, 0))
cov_50 <- interval_coverage(observed, predicted, quantile_level, interval_range = c(50))
cov_80 <- interval_coverage(observed, predicted, quantile_level, interval_range = c(80))
manual <- 0.5 * (cov_50 - 0.5) + 0.5 * (cov_80 - 0.8)
expect_equal(
interval_coverage_deviation(observed, predicted, quantile_level),
manual
)
expect_warning(
interval_coverage_deviation(
observed, predicted, c(quantile_level[-4], 0.76)
),
"all quantiles must form central symmetric prediction intervals"
)
})
# ============================================================================ #
# `bias_quantile` ============================================================ #
# ============================================================================ #
test_that("bias_quantile() works as expected", {
predicted <- c(1, 2, 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
0
)
predicted <- c(0, 1, 2)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
-0.8
)
predicted <- c(
705.500, 1127.000, 4006.250, 4341.500, 4709.000, 4821.996,
5340.500, 5451.000, 5703.500, 6087.014, 6329.500, 6341.000,
6352.500, 6594.986, 6978.500, 7231.000, 7341.500, 7860.004,
7973.000, 8340.500, 8675.750, 11555.000, 11976.500
)
quantile_level <- c(0.01, 0.025, seq(0.05, 0.95, 0.05), 0.975, 0.99)
observed <- 8062
expect_equal(bias_quantile(observed, predicted, quantile_level), -0.8)
})
test_that("bias_quantile handles matrix input", {
observed <- seq(10, 0, length.out = 4)
predicted <- matrix(1:12, ncol = 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed, predicted, quantiles),
c(-1.0, -0.8, 0.8, 1.0)
)
})
test_that("bias_quantile() handles vector that is too long", {
predicted <- c(NA, 1, 2, 3)
quantiles <- c(0.1, 0.5, 0.9)
expect_error(
bias_quantile(observed = 2, predicted, quantiles),
"Assertion on 'quantile_level' failed: Must have length 4, but has length 3."
)
})
test_that("bias_quantile() handles NA values", {
predicted <- c(NA, 1, 2)
quantiles <- c(0.1, 0.5, 0.9)
expect_equal(
bias_quantile(observed = 2, predicted, quantiles),
-0.8
)
predicted <- c(0, 1, 2)
quantiles <- c(0.1, 0.5, NA)
expect_equal(
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", {
expect_error(
bias_quantile(observed = 2, numeric(0), numeric(0)),
"Assertion on 'quantile_level' failed: Must have length >= 1, but has length 0"
)
})
test_that("bias_quantile() returns correct bias if value below the median", {
predicted <- c(1, 2, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 1, predicted, quantiles), 0.8)
)
})
test_that("bias_quantile() returns correct bias if value above median", {
predicted <- c(1, 2, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 5, predicted, quantiles), -0.8)
)
})
test_that("bias_quantile() returns correct bias if value at the median", {
predicted <- c(1, 2, 3, 4)
quantiles <- c(0.1, 0.3, 0.5, 0.7)
expect_equal(bias_quantile(observed = 3, predicted, quantiles), 0)
})
test_that("bias_quantile() returns 1 if true value below min prediction", {
predicted <- c(2, 3, 4, 5)
quantiles <- c(0.1, 0.3, 0.7, 0.9)
suppressMessages(
expect_equal(bias_quantile(observed = 1, predicted, quantiles), 1)
)
})
test_that("bias_quantile() returns -1 if true value above max prediction", {
predicted <- c(1, 2, 3, 4)
quantiles <- c(0.1, 0.3, 0.5, 0.7)
expect_equal(bias_quantile(observed = 6, predicted, quantiles), -1)
})
test_that("bias_quantile(): quantiles must be between 0 and 1", {
predicted <- 1:4
# Failing example
quantiles <- c(-0.1, 0.3, 0.5, 0.8)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Assertion on 'quantile_level' failed: Element 1 is not >= 0."
)
# Passing counter example
quantiles <- c(0.1, 0.3, 0.5, 0.8)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile(): quantiles must be increasing", {
predicted <- 1:4
# Failing example
quantiles <- c(0.8, 0.3, 0.5, 0.9)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Predictions must not be decreasing with increasing quantile level"
)
# Passing counter example
quantiles <- c(0.3, 0.5, 0.8, 0.9)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile(): predictions must be increasing", {
predicted <- c(1, 2, 4, 3)
quantiles <- c(0.1, 0.3, 0.5, 0.9)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Predictions must not be decreasing with increasing quantile level"
)
expect_silent(bias_quantile( observed = 3, 1:4, quantiles))
})
test_that("bias_quantile(): quantile levels must be unique", {
predicted <- 1:4
# Failing example
quantiles <- c(0.3, 0.3, 0.5, 0.8)
expect_error(
bias_quantile(observed = 3, predicted, quantiles),
"Assertion on 'quantile_level' failed: Contains duplicated values, position 2."
)
# Passing example
quantiles <- c(0.3, 0.5, 0.8, 0.9)
expect_silent(bias_quantile(observed = 3, predicted, quantiles))
})
test_that("bias_quantile only produces one message", {
expect_message(
bias_quantile(observed, predicted[, -3], quantile_level[-3]),
"Median not available, computing bias as mean of the two innermost quantiles in order to compute bias."
)
})
test_that("bias_quantile() works with point forecasts", {
predicted <- 1
observed <- 1
quantile_level <- 0.5
expect_equal(bias_quantile(observed, predicted, quantile_level), 0)
})
# `interpolate_median` ======================================================= #
test_that("interpolation in `interpolate_median` works", {
predicted <- c(1, 10)
observed <- 15
quantile_level <- c(0.4, 0.6)
# median is missing, symmetric quantile levels given
expect_equal(interpolate_median(predicted, quantile_level), mean(predicted))
quantile_level <- c(0.1, 0.9)
expect_equal(interpolate_median(predicted, quantile_level), mean(predicted))
1 + 0.4/0.8 * (10 - 1)
# asymmetric quantile levels given
quantile_level <- c(0.3, 0.6)
expect_equal(interpolate_median(predicted, quantile_level), 1 + 0.2/0.3 * (10 - 1))
quantile_level <- c(0.2, 0.6)
expect_equal(interpolate_median(predicted, quantile_level), 1 + 0.3/0.4 * (10 - 1))
})
# ============================================================================ #
# `ae_median_quantile` ======================================================= #
# ============================================================================ #
test_that("ae_median_quantile() works as_expected", {
observed <- 1:30
predicted_values <- matrix(2:31)
expect_equal(
ae_median_quantile(observed, predicted_values, quantile_level = 0.5),
as.numeric(predicted_values - observed)
)
# using a vector as input for predicted does not work if there are several
# observed values
expect_error(
ae_median_quantile(observed, as.numeric(predicted_values), quantile_level = 0.5),
"Assertion on 'predicted' failed: Must be of type 'matrix', not 'double'."
)
# it does work if there is only a single observed value - in this case
# the vector gets treated as one forecast - which means that we need to
# supply multiple quantile levels now, because it assumes that single forecast
# is represented by several quantiles
expect_equal(
ae_median_quantile(observed[1], as.numeric(predicted_values)[1:3], quantile_level = c(0.1, 0.5, 0.7)),
2
)
# 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.'
)
dim(1:10)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.