scoring (version 0.4)

calcscore: Calculate Scores Under A Specific Rule

Description

Given parameters of a scoring rule family, calculate scores for probabilistic forecasts and associated outcomes.

Usage

## S3 method for class 'formula':
calcscore(object, fam="pow", param, data, bounds=NULL,
                  reverse=FALSE, ...)
## S3 method for class 'default':
calcscore(object, outcome, fam="pow",
                  param=c(2,rep(1/max(2,NCOL(forecast)),max(2,NCOL(forecast)))),
                  bounds=NULL, reverse=FALSE, ...)

Arguments

object

an object of class "formula", of the form outcome ~ forecast (see details). Alternatively, a matrix of forecasts, with observations in rows and forecast alternatives in columns. For two-alternative forecasts, this

outcome

a vector of outcomes. For each row of the forecast matrix, outcome should contain an entry reflecting the column number associated with the event that occurred.

fam

scoring rule family. pow (default) is the power family, beta is the beta family, sph is the pseudospherical family.

param

for family beta, a numeric vector of length 2 containing the scoring rule family parameters. For other families, a numeric vector of length NCOL(forecast)+1, containing first the family parameter gamma and then th

data

an optional data frame or list containing the variables in the model. If not found in data, the variables are taken from the environment from which calcscore is called.

bounds

a vector of length 2 corresponding to the desired minimum value and maximum value of the scoring rule, respectively. Entries of NA imply that the minimum and/or maximum bound will not be modified from the natural, fam

reverse

if FALSE (default), smaller scores imply better forecasts. If TRUE, larger scores imply better forecasts.

...

Additional arguments.

Value

• calcscore returns a numeric vector that has length equal to length(outcome), containing scores under the selected scoring rule.

Details

The formula is of the form outcome ~ forecast, where forecast describes the column(s) containing forecasts associated with the possible outcomes. Multiple columns are separated by +. outcome is always a vector describing the outcome associated with each forecast. These should be coded 1, 2, ..., reflecting the column associated with the outcome (see examples).

For events with only two alternatives, one can take a shortcut and supply only forecasts and baseline parameters (for families pow and sph) associated with a single outcome. In this case, the outcome vector should contain zeros and ones, where `one' means that the forecasted alternative occurred.

When the bounds argument is supplied, the code attempts to scale the scores so that the maximum score is bounds[2] and the minimum score is bounds[1]. This scaling cannot be accomplished when the scoring rule allows scores of infinity (the log score is the most common case here). If reverse=TRUE, the bounds are applied after the reversal (so that the supplied lower bound reflects the worst score and upper bound reflects the best score).

References

Buja, A., Stuetzle, W., & Shen, Y. (2005). Loss functions for binary class probability estimation and classification: Structure and applications. (Obtained from http://stat.wharton.upenn.edu/~buja/PAPERS/)

Jose, V. R. R., Nau, R. F., & Winkler, R. L. (2008). Scoring rules, generalized entropy, and utility maximization. Operations Research, 56, 1146--1157.

Jose, V. R. R., Nau, R. F., & Winkler, R. L. (2009). Sensitivity to distance and baseline distributions in forecast evaluation. Management Science, 55, 582--590.

Merkle, E. C. & Steyvers, M. (in press). Choosing a strictly proper scoring rule. Decision Analysis.

See Also

plotscore

Examples

## Brier scores for two alternatives, with bounds of 0 and 1
data("WorldEvents")
scores <- calcscore(answer ~ forecast, fam="beta",
                    param=c(1,1), data=WorldEvents,
                    bounds=c(0,1))

## Calculate Brier scores manually
scores.man <- with(WorldEvents, (forecast - answer)^2)

## Comparison
all.equal(scores, scores.man)

## Average Brier score for each forecaster
with(WorldEvents, tapply(scores, forecaster, mean))

## Brier scores for 3 alternatives, with bounds of 0 and 1
data("WeatherProbs")
scores2 <- calcscore(tcat ~ tblw + tnrm + tabv, fam="pow",
                     param=c(2,rep(1/3,3)), data=WeatherProbs,
                     bounds=c(0,1))

## Spherical scores for 3 alternatives, reversed so 0 is worst and
## 1 is best
scores3 <- calcscore(tcat ~ tblw + tnrm + tabv, fam="sph",
                     param=c(2,rep(1/3,3)), data=WeatherProbs,
                     bounds=c(0,1), reverse=TRUE)