Negative values of CRPS?

[mbojan]

I'm getting the following when trying to score predictions from a negative binomial GLM:

dat <- list(
  observed = c(
    65L, 62L, 48L, 47L, 60L, 75L, 62L, 50L, 49L, 
    54L, 62L, 51L, 57L, 56L, 48L, 47L, 43L, 35L, 53L, 34L, 58L, 38L, 
    56L, 42L, 58L, 55L, 51L, 54L, 64L, 51L, 45L, 59L, 45L, 46L, 67L, 
    44L, 51L, 47L, 50L, 56L, 69L, 56L, 35L, 39L, 50L, 45L, 43L, 51L, 
    47L, 53L, 36L, 30L
  ), 
  pred_means = c(
    48.3099006480877, 48.1856853530793, 
    48.0249249252662, 47.8581559836325, 47.7155235009397, 47.6268375218476, 
    47.6218179266744, 47.724978436259, 47.9185569151005, 48.1648674172695, 
    48.4256057306771, 48.6616929438018, 48.833317728955, 48.9007310037522, 
    48.8620539654412, 48.7821309824475, 48.7324016131518, 48.7840863337008, 
    49.0091496602714, 49.4823260083412, 50.2675286039073, 51.3049811917389, 
    52.4702988591609, 53.6283020516166, 54.6318059743995, 55.3238719063732, 
    55.5453815997929, 55.2375969406439, 54.5183299747122, 53.5268537918377, 
    52.3989708593686, 51.2620683576382, 50.2333662682769, 49.4087504392848, 
    48.7853128755403, 48.3188875640644, 47.9682808060746, 47.6944186604847, 
    47.4595844049793, 47.2271670054595, 46.9786501957888, 46.7263257219719, 
    46.4853115157984, 46.2704191049954, 46.0962135356727, 45.9771198543237, 
    45.9245710565769, 45.9269448553235, 45.9619878567896, 46.0074298534593, 
    46.0409403503701, 46.0401323330852
  ), 
  size = 249.137522940415
)

scoringRules::crps(
  dat$observed,
  family = "negative-binomial",
  mu = dat$pred_means,
  size = dat$size
)
#> Warning in crps.numeric(dat$observed, family = "negative-binomial", mu =
#> dat$pred_means, : Negative CRPS values. Check parameter combinations and contact
#> package maintainer(s).
#>  [1] -7.128161e+23  5.192899e+23  6.433863e+23 -1.274805e+23 -6.621426e+22
#>  [6]  8.679541e+23 -1.282874e+23  3.759245e+23  1.074203e+23  1.297433e+24
#> [11]  2.294787e+23 -4.796344e+23  7.317828e+23 -2.552954e+23 -3.320416e+22
#> [16]  9.255501e+22  9.219190e+23  4.415978e+22  1.329747e+24  2.842438e+24
#> [21]  3.246396e+24  6.936429e+24  2.548535e+00  7.520443e+00  2.523148e+00
#> [26]  1.915409e+00  2.806529e+00  1.961657e+00  5.981568e+00  2.137097e+00
#> [31]  4.305351e+00  1.150249e+25  4.810089e+24  2.194627e+24  6.986858e+22
#> [36] -8.343210e+23 -2.429890e+23  8.616790e+22  8.886526e+22 -5.083535e+22
#> [41]  1.869346e+23  4.855898e+22  1.920124e+22 -1.208002e+23  5.227473e+21
#> [46]  4.265401e+22 -2.974556e+22 -2.418369e+22  5.379786e+22 -3.412752e+21
#> [51] -2.363378e+22  8.044060e+22

[aijordan]

Oh, that's not great. Thanks for the report.

I'll make a suggestion in the coming days, when I find the time: The problem is the automatic dispatch of an appropriate expression for the Gaussian hypergeometric function in the hypergeo package, which never claimed to guarantee accurate results for large parameter values (e.g., a size parameter of approx 250). Obviously, these values are reasonable in negative binomial GLMs, and we need to address this.

[mbojan]

Thanks @aijordan . I'll appreciate the suggestion. Let me know if you need any additional info.

[mbojan]

@aijordan can you point me to the paper that inspired your implementation? I went through Gneiting, T. and A.E. Raftery (2007) and Gneiting, T. and M. Katzfuss (2014) that are cited on the man page, but I don't understand the step you made to use the hypergeometric function.

[mbojan]

I decided, for now, to get away with scoringRules::crps_sample() and

dat <- list(
  observed = c(
    65L, 62L, 48L, 47L, 60L, 75L, 62L, 50L, 49L, 
    54L, 62L, 51L, 57L, 56L, 48L, 47L, 43L, 35L, 53L, 34L, 58L, 38L, 
    56L, 42L, 58L, 55L, 51L, 54L, 64L, 51L, 45L, 59L, 45L, 46L, 67L, 
    44L, 51L, 47L, 50L, 56L, 69L, 56L, 35L, 39L, 50L, 45L, 43L, 51L, 
    47L, 53L, 36L, 30L
  ), 
  pred_means = c(
    48.3099006480877, 48.1856853530793, 
    48.0249249252662, 47.8581559836325, 47.7155235009397, 47.6268375218476, 
    47.6218179266744, 47.724978436259, 47.9185569151005, 48.1648674172695, 
    48.4256057306771, 48.6616929438018, 48.833317728955, 48.9007310037522, 
    48.8620539654412, 48.7821309824475, 48.7324016131518, 48.7840863337008, 
    49.0091496602714, 49.4823260083412, 50.2675286039073, 51.3049811917389, 
    52.4702988591609, 53.6283020516166, 54.6318059743995, 55.3238719063732, 
    55.5453815997929, 55.2375969406439, 54.5183299747122, 53.5268537918377, 
    52.3989708593686, 51.2620683576382, 50.2333662682769, 49.4087504392848, 
    48.7853128755403, 48.3188875640644, 47.9682808060746, 47.6944186604847, 
    47.4595844049793, 47.2271670054595, 46.9786501957888, 46.7263257219719, 
    46.4853115157984, 46.2704191049954, 46.0962135356727, 45.9771198543237, 
    45.9245710565769, 45.9269448553235, 45.9619878567896, 46.0074298534593, 
    46.0409403503701, 46.0401323330852
  ), 
  size = 249.137522940415
)

scoringRules::crps(
  dat$observed,
  family = "negative-binomial",
  mu = dat$pred_means,
  size = dat$size
)
#> Warning in crps.numeric(dat$observed, family = "negative-binomial", mu =
#> dat$pred_means, : Negative CRPS values. Check parameter combinations and contact
#> package maintainer(s).
#>  [1] -7.128161e+23  5.192899e+23  6.433863e+23 -1.274805e+23 -6.621426e+22
#>  [6]  8.679541e+23 -1.282874e+23  3.759245e+23  1.074203e+23  1.297433e+24
#> [11]  2.294787e+23 -4.796344e+23  7.317828e+23 -2.552954e+23 -3.320416e+22
#> [16]  9.255501e+22  9.219190e+23  4.415978e+22  1.329747e+24  2.842438e+24
#> [21]  3.246396e+24  6.936429e+24  2.548535e+00  7.520443e+00  2.523148e+00
#> [26]  1.915409e+00  2.806529e+00  1.961657e+00  5.981568e+00  2.137097e+00
#> [31]  4.305351e+00  1.150249e+25  4.810089e+24  2.194627e+24  6.986858e+22
#> [36] -8.343210e+23 -2.429890e+23  8.616790e+22  8.886526e+22 -5.083535e+22
#> [41]  1.869346e+23  4.855898e+22  1.920124e+22 -1.208002e+23  5.227473e+21
#> [46]  4.265401e+22 -2.974556e+22 -2.418369e+22  5.379786e+22 -3.412752e+21
#> [51] -2.363378e+22  8.044060e+22

z <- sapply(
  seq(along=dat$observed),
  function(i) {
    scoringRules::crps_sample(
      y = dat$observed[i],
      dat = rnbinom(10000, size=dat$size, mu = dat$pred_means[i])
    )
  }
)
z
#>  [1] 12.495670  9.723042  1.760610  1.770754  8.536062 23.141016 10.300841
#>  [8]  2.115734  1.864141  3.479566  9.527966  2.152530  5.157453  4.241641
#> [15]  1.790685  1.917909  3.213451  9.754882  2.730997 11.241431  4.641654
#> [22]  9.254674  2.554329  7.460273  2.487824  1.917368  2.888156  1.962453
#> [29]  5.919695  2.118289  4.295258  4.725922  3.058588  2.267781 14.017569
#> [36]  2.624592  2.275088  1.770857  2.105888  5.556528 17.920134  5.951265
#> [43]  7.534502  4.244500  2.640825  1.744088  2.124471  3.170937  1.798122
#> [50]  4.266090  6.340067 11.827803

What do you think @aijordan ? Are these viable approximations? Do you have other suggestions?

[aijordan]

@mbojan, let me suggest this slight variation of your solution for now, without RNG:

sapply(
  seq(along=dat$observed),
  function(i) {
    scoringRules::crps_sample(
      y = dat$observed[i],
      dat = 0:150,
      w = dnbinom(0:150, size=dat$size, mu = dat$pred_means[i])
    )
  }
)
#>   [1] 12.524908  9.811189  1.764168  1.779142  8.441514 23.132536 10.355357
#>   [8]  2.083245  1.850250  3.582602  9.582249  2.116417  5.081308  4.352307
#> [15]  1.800542  1.903845  3.316265  9.631608  2.681845 11.223429  4.766892
#> [22]  9.090524  2.548535  7.520443  2.523148  1.915409  2.806529  1.961657
#> [29]  5.981568  2.137097  4.305351  4.764429  3.066962  2.312115 13.983689
#> [36]  2.635610  2.310317  1.766204  2.150713  5.545419 17.828083  5.938077
#> [43]  7.604773  4.256486  2.611230  1.748063  2.109028  3.158801  1.803828
#> [50]  4.294543  6.368936 11.923287

Just make sure that the support is large enough, which is easily the case with 0:150 for your example:

print(pnbinom(150, size = dat$size, mu = dat$pred_means, lower.tail = FALSE), digits = 1)
#>   [1] 3e-25 2e-25 2e-25 1e-25 1e-25 8e-26 8e-26 1e-25 1e-25 2e-25 3e-25 5e-25
#> [13] 7e-25 8e-25 8e-25 7e-25 6e-25 7e-25 1e-24 2e-24 8e-24 5e-23 3e-22 2e-21
#> [25] 8e-21 2e-20 3e-20 2e-20 7e-21 2e-21 3e-22 4e-23 8e-24 2e-24 7e-25 3e-25
#> [37] 2e-25 9e-26 6e-26 4e-26 2e-26 2e-26 1e-26 6e-27 5e-27 4e-27 3e-27 3e-27
#> [49] 4e-27 4e-27 4e-27 4e-27

[mbojan]

Thanks @aijordan !

[FK83]

Thanks, @mbojan and @aijordan! Just for completeness, the formulas for the CRPS are listed in the appendix of the paper accompanying the package (Jordan et al., Journal of Statistical Software 90, 2019, https://doi.org/10.18637/jss.v090.i12).

[mbojan]

Thanks for the reference @FK83 . Are you sure you want to close this? The issue with crps() still stands.

[FK83]

Agreed, that was too fast - just reopened the issue.