R/calibration.R

In surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

################################################################################
### Calibration tests for count data based on proper scoring rules
### Reference: Wei and Held (2014), Test, 23, 787-805
###
### Copyright (C) 2015,2018 Sebastian Meyer
###
### This file is part of the R package "surveillance",
### free software under the terms of the GNU General Public License, version 2,
### a copy of which is available at https://www.R-project.org/Licenses/.
################################################################################


## perform a calibration test given observations x
## with Poisson (size = NULL) or NegBin predictions

calibrationTest.default <- function (x, mu, size = NULL,
                                     which = c("dss", "logs", "rps"),
                                     tolerance = 1e-4, method = 2, ...)
{
    stopifnot(x >= 0, mu > 0, is.null(size) || all(size > 0))

    ## calculate scores
    which <- match.arg(which)
    score <- do.call(which, args = alist(x = x, mu = mu, size = size))

    ## calculate z-statistic
    z <- calibrationZ(score, mu, size, which, tolerance, method)

    ## calculate two-sided p-value
    p <- 2 * pnorm(-abs(z))

    ## construct an object of class "htest"
    res <- list(
        method = paste0("Calibration Test for Count Data (based on ",
            toupper(which), ")"),
        data.name = deparse(substitute(x)),
        statistic = c("z" = z),
        parameter = c("n" = length(x)),
        p.value = p
    )
    class(res) <- "htest"
    res
}

## compute the calibration z-statistic given the computed scores
calibrationZ <- function (score, mu, size = NULL,
                          which = c("dss", "logs", "rps"),
                          tolerance = 1e-4, method = 2)
{
    stopifnot(method %in% 1:2)
    ## expectation and variance of score for given predictive distribution
    EV <- score_EV(mu, size, tolerance, which)
    ## calculate the z-statistic
    z <- do.call(paste0("zScore", method),
                 args = alist(score, EV[[1L]], EV[[2L]]))
    z
}

## compute the calibration z-statistic and p-value
## from a set of scores and their null expectations and variances
zScore1 <- function (score, E0, V0)
{
    n <- length(score)

    ## emean <- mean(E0)
    ## varmean <- sum(V0) / n^2
    ## (mean(score) - emean) / sqrt(varmean)

    sum(score - E0) / sqrt(sum(V0))
}

## alternative z-statistic Z*
zScore2 <- function (score, E0, V0)
{
    n <- length(score)
    sum((score - E0) / sqrt(V0)) / sqrt(n)
}