Subsetting reliability diagram objects

Description

Subsetting reliability diagram objects

Usage

## S3 method for class 'reliabilitydiag'
x[i]

Arguments

Value

an object inheriting from the class 'reliabilitydiag'.

See Also

c.reliabilitydiag.

Examples

data("precip_Niamey_2016", package = "reliabilitydiag")

r <- reliabilitydiag(
  precip_Niamey_2016[c("Logistic", "EMOS")],
  y = precip_Niamey_2016$obs
)
length(r)
r[1]
r["EMOS"]

Coerce to a reliability diagram

Description

Coerce numeric vectors, data frames, or anything else that can be coerced by as.data.frame to a data frame of prediction values, into an object inheriting from the 'reliabilitydiag' class.

Usage

as.reliabilitydiag(x, ...)

is.reliabilitydiag(x)

## S3 method for class 'reliabilitydiag'
as.reliabilitydiag(x, y = NULL, r = NULL, tol = sqrt(.Machine$double.eps), ...)

## Default S3 method:
as.reliabilitydiag(
  x,
  y = NULL,
  r = NULL,
  xtype = NULL,
  xvalues = NULL,
  .name_repair = "unique",
  region.level = 0.9,
  region.method = NULL,
  region.position = "diagonal",
  n.boot = 100,
  ...
)

## S3 method for class 'data.frame'
as.reliabilitydiag(
  x,
  y = NULL,
  r = NULL,
  xtype = NULL,
  xvalues = NULL,
  .name_repair = "unique",
  region.level = 0.9,
  region.method = NULL,
  region.position = "diagonal",
  n.boot = 100,
  ...
)

Arguments

Value

as.reliabilitydiag returns a 'reliabilitydiag' object.

is.reliabilitydiag returns TRUE if its argument is a reliability diagram, that is, has "reliabilitydiag" among its classes, and FALSE otherwise.

See Also

reliabilitydiag

Combining reliability diagram objects

Description

Combine two or more 'reliabilitydiag' objects that are based on the same observations. Other objects are coerced by as.reliabilitydiag before combination.

Usage

## S3 method for class 'reliabilitydiag'
c(
  ...,
  tol = sqrt(.Machine$double.eps),
  xtype = NULL,
  xvalues = NULL,
  region.level = 0.9,
  region.method = NULL,
  region.position = "diagonal",
  n.boot = 100
)

Arguments

Value

an object inheriting from the class 'reliabilitydiag'.

See Also

as.reliabilitydiag, [.reliabilitydiag.

Examples

data("precip_Niamey_2016", package = "reliabilitydiag")

X <- precip_Niamey_2016[c("EMOS", "ENS")]
Y <- precip_Niamey_2016$obs
r0 <- reliabilitydiag0(Y)
r1 <- c(r0, X, EPC = precip_Niamey_2016$EPC, region.level = NA)
r1
c(r1, reliabilitydiag(Logistic = precip_Niamey_2016$Logistic, y = Y))

Miscalibration Test

Description

(experimental)

Usage

miscalibration_test(x, ...)

## S3 method for class 'reliabilitydiag'
miscalibration_test(x, ...)

## S3 method for class 'numeric'
miscalibration_test(x, y, ...)

Arguments

Value

returns a 'tibble' with entries

Plotting reliability diagram objects

Description

Using the ggplot2 package to visually diagnose the reliability of prediction methods that issue probability forecasts.

Usage

## S3 method for class 'reliabilitydiag'
plot(x, ...)

## S3 method for class 'reliabilitydiag'
autoplot(
  object,
  ...,
  type = c("miscalibration", "discrimination"),
  colour = "red",
  params_histogram = NULL,
  params_ggMarginal = NULL,
  params_ribbon = NULL,
  params_diagonal = NULL,
  params_vsegment = NULL,
  params_hsegment = NULL,
  params_CEPline = NULL,
  params_CEPsegment = NULL,
  params_CEPpoint = NULL
)

## S3 method for class 'reliabilitydiag'
autolayer(
  object,
  ...,
  type = c("miscalibration", "discrimination"),
  colour = "red",
  params_histogram = NA,
  params_ggMarginal = NA,
  params_ribbon = NA,
  params_diagonal = NA,
  params_vsegment = NA,
  params_hsegment = NA,
  params_CEPline = NA,
  params_CEPsegment = NA,
  params_CEPpoint = NA
)

Arguments

Details

plot always sends a plot to a graphics device, wheres autoplot behaves as any ggplot() + layer() combination. That means, customized plots should be created using autoplot and autolayer.

Three sets of default parameter values are used:

• If multiple predictions methods are compared, then only the most necessary information to determine reliability are displayed.

• For a single prediction method and type = "miscalibration", the focus lies on the deviation from the diagonal including uncertainty quantification.

• For a single prediction method and type = "discrimination", the focus lies on the PAV transformation and the resulting marginal distribution. A concentration of CEP values near 0 or 1 suggest a high potential predictive ability of a prediction method.

Setting any of the params_* arguments to NA disables that layer.

Default parameter values if length(object) > 1, where the internal variable forecast is used as grouping variable:

Default parameter values for type = "miscalibration" if length(object) == 1:

Default parameter values for type = "discrimination" if length(object) == 1:

Value

An object inheriting from class 'ggplot'.

Examples

data("precip_Niamey_2016", package = "reliabilitydiag")
r <- reliabilitydiag(
  precip_Niamey_2016[c("Logistic", "EMOS", "ENS", "EPC")],
  y = precip_Niamey_2016$obs,
  region.level = NA
)

# simple plotting
plot(r)

# faceting using the internal grouping variable 'forecast'
autoplot(r, params_histogram = list(colour = "black", fill = NA)) +
  ggplot2::facet_wrap("forecast")

# custom color scale for multiple prediction methods
cols <- c(Logistic = "red", EMOS = "blue", ENS = "darkgreen", EPC = "orange")
autoplot(r) +
  ggplot2::scale_color_manual(values = cols)

# default reliability diagram type with a title
rr <- reliabilitydiag(
  EMOS = precip_Niamey_2016$EMOS,
  r = r,
  region.level = 0.9
)
autoplot(rr) +
  ggplot2::ggtitle("Reliability diagram for EMOS method")

# using defaults for discrimination diagrams
p <- autoplot(r["EMOS"], type = "discrimination")
print(p, newpage = TRUE)

# ggMarginal needs to be called after adding all custom layers
p <- autoplot(r["EMOS"], type = "discrimination", params_ggMarginal = NA) +
  ggplot2::ggtitle("Discrimination diagram for EMOS method")
p <- ggExtra::ggMarginal(p, type = "histogram")
print(p, newpage = TRUE)

# the order of the layers can be changed
autoplot(rr, colour = "black", params_ribbon = NA) +
  autolayer(rr, params_ribbon = list(fill = "red", alpha = .5))

Precipitation forecasts and observations at Niamey, Niger in July to September 2016

Description

A data set containing 24-hour ahead daily probability of precipitation forecasts of four forecasting methods and corresponding observations of precipitation occurrence.

For a detailed description of the four prediction methods, see Vogel et al (2021).

Usage

precip_Niamey_2016

Format

A data frame with 92 rows and 6 variables:

date

a date from "2016-07-01" to "2016-09-30" in Date format.

Logistic

prediction based on logistic regression, as a probability.

EMOS

prediction based on EMOS method, as a probability.

ENS

prediction based on ECMWF raw ensemble, as a probability.

EPC

prediction based on EPC method, as a probability.

obs

observation, indicator variable where 1 represents the occurrence of precipitation.

Source

Vogel P, Knippertz P, Gneiting T, Fink AH, Klar M, Schlueter A (2021). "Statistical forecasts for the occurrence of precipitation outperform global models over northern tropical Africa." Geophysical Research Letters, 48, e2020GL091022. doi:10.1029/2020GL091022.

This data set contains modified historic products from the European Center for Medium-Range Weather Forecasts (ECMWF, https://www.ecmwf.int/), specifically: ensemble forecasts of precipitation that have been summarized to a probability of precipitation (column ENS), and historical observations for the occurence of precipitation (column obs). The ECMWF licenses the use of expired real-time data products under the Creative Commons Attribution 4.0 International (CC BY 4.0, https://creativecommons.org/licenses/by/4.0/).

Printing reliability diagram objects

Description

Printing methods for 'reliabilitydiag' and 'summary.reliabilitydiag' objects.

Usage

## S3 method for class 'reliabilitydiag'
print(x, ...)

## S3 method for class 'summary.reliabilitydiag'
print(x, ...)

Arguments

Details

print.reliabilitydiag always sends a plot to the current graphics device and prints a summary to the console.

print.summary.reliabilitydiag prints the summary output to the console.

Value

Invisibly returns x.

See Also

autoplot.reliabilitydiag, summary.reliabilitydiag

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

ggplot2

autolayer, autoplot

Reliability diagram object

Description

Documentation of the 'reliabilitydiag' object, and its constructors.

Usage

reliabilitydiag(
  ...,
  y = NULL,
  r = NULL,
  tol = sqrt(.Machine$double.eps),
  xtype = NULL,
  xvalues = NULL,
  region.level = 0.9,
  region.method = NULL,
  region.position = "diagonal",
  n.boot = 100
)

reliabilitydiag0(y)

Arguments

Details

reliabilitydiag constructs and returns an object inheriting from the class 'reliabilitydiag'. Each object passed via ... is coerced by the methods described in as.reliabilitydiag, and then concatenated by c.reliabilitydiag.

reliabilitydiag0 constructs an empty 'reliabilitydiag' object from the response values.

If any of the arguments region.level, region.method, or region.position is NA, then the uncertainty quantification in terms of consistency/confidence regions is skipped.

Consistency regions are determined under the assumption of calibration of the original predictions, that is, perfectly reliable forecasts such that P(Y = 1|X) = X. Consistency regions are therefore positioned around values on the diagonal (set region.position to "diagonal").

For confidence regions, calibration is enforced by using the PAV-recalibrated predictions for uncertainty quantification, that is, it is assumed that P(Y = 1|X) = PAV(X). Confidence regions are therefore positioned around the estimated conditional exceedence probability (CEP) line (set region.position to "estimate").

When region.method is "resampling", then the original forecast-observations pairs are bootstrapped n.boot times. For each bootstrap sample, new observations are drawn under the respective assumption (consistency or confidence). Then PAV-recalibration with those new observations is performed on each bootstrap sample, and pointwise lower and upper bounds are calculated across the resulting CEP lines.

When region.method is "discrete_asymptotics" and region.position is "diagonal", a Gaussian approximation is used assuming \sqrt{n} * (EST(x) - x) has variance x(1-x), where x is an original prediction value, n is the observed number of predictions with value x, and EST(x) is the estimated CEP value at x.

When region.method is "continuous_asymptotics" and region.position is "diagonal", a Chernoff approximation is used for (n * f(x) / (4 * x * (1- x)))^{(1/3)} * (EST(x) - x), where x is an original prediction value, n is the total number of observations, EST(x) is the estimated CEP value at x, and f(x) is the estimated value of the density of the original prediction values. This density is estimated using the bde package: We use Chen's beta kernel density estimator (see bde).

Value

reliabilitydiag returns a 'reliabilitydiag' object, which is a named list-type vector class with the attribute y containing the values supplied to the input argument y, that is, the numeric vector of response values to be predicted. The length is given by the number of prediction methods detected from the supplied objects.

reliabilitydiag0 returns an empty 'reliabilitydiag' object with attribute y.

Each entry of a 'reliabilitydiag' object (corresponding to a single prediction method) is itself a list with the following entries

Each cases tibble comprises the forecast-observation pairs of the given prediction method. It is arranged in increasing order of x and has columns

Each bins tibble contains PAV-recalibration information, and has columns

Each regions tibble contains the uncertainty quantification information, and has columns

Each xinfo list has entries

See Also

c.reliabilitydiag, [.reliabilitydiag, plot.reliabilitydiag.

See summary.reliabilitydiag for a decomposition of predictive performance into miscalibration, discrimination, and uncertainty.

Examples

data("precip_Niamey_2016", package = "reliabilitydiag")

# standard use with a data.frame
r <- reliabilitydiag(precip_Niamey_2016["EMOS"], y = precip_Niamey_2016$obs)
r

# no consistency/confidence regions
X <- precip_Niamey_2016$EMOS
Y <- precip_Niamey_2016$obs
r1 <- reliabilitydiag(X = X, y = Y, region.level = NA)
r1

# specify predictions via existing reliabilitydiag
r0 <- reliabilitydiag0(Y)
identical(r1, reliabilitydiag(X = X, r = r0, region.level = NA))

# only observation information is used from existing reliabilitydiag
X2 <- precip_Niamey_2016$ENS
r2 <- reliabilitydiag(X2 = X2, r = r, region.level = NA)
r3 <- reliabilitydiag(X2 = X2, r = r0, region.level = NA)
identical(r2, r3)

Decomposing scores into miscalibration, discrimination and uncertainty

Description

An object of class reliabilitydiag contains the observations, the original forecasts, and recalibrated forecasts given by isotonic regression. The function summary.reliabilitydiag calculates quantitative measures of predictive performance, miscalibration, discrimination, and uncertainty, for each of the prediction methods in relation to their recalibrated version.

Usage

## S3 method for class 'reliabilitydiag'
summary(object, ..., score = "brier")

Arguments

Details

Predictive performance is measured by the mean score of the original forecast values, denoted by S.

Uncertainty, denoted by UNC, is the mean score of a constant prediction at the value of the average observation. It is the highest possible mean score of a calibrated prediction method.

Discrimination, denoted by DSC, is UNC minus the mean score of the PAV-recalibrated forecast values. A small value indicates a low information content (low signal) in the original forecast values.

Miscalibration, denoted by MCB, is S minus the mean score of the PAV-recalibrated forecast values. A high value indicates that predictive performance of the prediction method can be improved by recalibration.

These measures are related by the following equation,

S = MCB - DSC + UNC.

Score decompositions of this type have been studied extensively, but the optimality of the PAV solution ensures that MCB is nonnegative, regardless of the chosen (admissible) scoring function. This is a unique property achieved by choosing PAV-recalibration.

If deviating from the Brier score as performance metric, make sure to choose a proper scoring rule for binary events, or equivalently, a scoring function with outcome space {0, 1} that is consistent for the expectation functional.

Value

A 'summary.reliability' object, which is also a tibble (see tibble::tibble()) with columns:

Examples

data("precip_Niamey_2016", package = "reliabilitydiag")
r <- reliabilitydiag(
  precip_Niamey_2016[c("Logistic", "EMOS", "ENS", "EPC")],
  y = precip_Niamey_2016$obs,
  region.level = NA
)
summary(r)
summary(r, score = function(y, x) (x - y)^2)