RAnEnExtra::verifyBrier

RAnEnExtra::verifyBrier computes the Brier score and its decomposition for each all lead times available. This function used the NCAR verification package.

verifyBrier(
  anen.ver,
  obs.ver,
  threshold,

    ensemble.func = stop("Ensemble function missing. If it is not an ensemble, set baseline = T"),
  baseline = F,
  ...
)

Arguments

anen.ver	A 4-dimensional array for analogs.
obs.ver	A 3-dimensional array for observations.
threshold	The numeric threshold for computing the brier score. Observation larger than or equal to the threshold will be converted to 1. Analog values will not be processed with this threshold value because it is assumed that the ensemble function ensemble.func will convert analog ensemble to probability values. Please see examples.
ensemble.func	A function that takes a vector as input and then converts the ensemble members (the 4th dimension of analogs) into a scalar. This scalar is usually a probability within [0, 1]. Please see examples.
baseline	TRUE if the anen.ver is an one-member forecast array. This is usually used for deterministic baseline model.
...	Extra parameters for the ensemble.func.

Examples

# Create synthetic datasets from different distributions
# with 10 stations, 2 test days, 5 lead times, and
# 50 ensemble members.
# 

anen.ver <- array(rnorm(5000), dim = c(10, 2, 5, 50))
obs.ver <- array(runif(100, min = -5, max = 5), dim = c(10, 2, 5))

# Set the threshold for observation. Values larger
# than or equal to this value will be converted to 1, otherwise 
# to 0.
#
threshold <- 0.5

# This is the function that takes the ensemble value vector 
# as input and output a probability. For example, here I count
# how many values in the member are larger than or equal to
# the split value. You can define your own function to compute
# a probability from the ensemble.
# 
ensemble.func <- function(v, split) {return(
  length(which(v >= split)) / length(v)
)}

# Calculate the Brier score.
# Don't forget to pass the additional threshold to your ensemble
# function. Otherwise you will receive an error complaining 
# about a missing argument.
# 
score <- verifyBrier(anen.ver, obs.ver, threshold, ensemble.func, split = threshold)

print(score)
#>          bs          ss reliability        resol
#> FLT1 0.3115 -0.29791667  0.12364286 0.0521428571
#> FLT2 0.3655 -0.60659341  0.13811905 0.0001190476
#> FLT3 0.2665 -0.06600000  0.06150000 0.0450000000
#> FLT4 0.2345  0.02291667  0.01783333 0.0233333333
#> FLT5 0.2955 -0.18200000  0.09073810 0.0452380952
#> All  0.2947 -0.18305901  0.05670971 0.0111097132

if (all(c('ggplot2', 'reshape2') %in% installed.packages())) {
  # Use ggplot2
  library(reshape2)
  library(ggplot2)

  # Unpivot the score table
  colnames(score) <- c('Brier', 'Skill', "Reliability", 'Resolution')
  melted <- melt(data = score)
  colnames(melted) <- c('FLT', 'Type', 'Value')

  # Generate ggplot
  ggplot(data = melted) +
    theme_bw() +
    geom_bar(mapping = aes(x = FLT, y = Value, group = Type, fill = Type),
             stat = 'identity', position = 'dodge') +
    scale_fill_brewer(palette = 'Paired', direction = -1) +
    geom_text(mapping = aes(x = FLT, y = Value, label = round(Value, digits = 3),
                            group = Type), position = position_dodge(0.9), vjust = -0.2) +
    labs(y = 'Score') +
    theme(legend.position = 'top')

} else {
  # Use the base graphics
  barplot(score[, 'bs'], ylab = 'Brier')
}