Empirical-Bayes shrinkage of hot-spot counts

The places with the highest counts in one period tend to fall in the next by chance alone: a count y = theta + e is a true rate plus noise of variance v (for Poisson counts, the mean count, the default). The linear shrinkage estimate (1 - B) y + B ybar toward the weighted mean has expected loss (1 - B)^2 S_e + B^2 S_theta when the noise has weighted mean 0 and is uncorrelated with the rates, which is smallest at B* = S_e / (S_e + S_theta), where it never exceeds the loss of the raw counts (James and Stein 1961; Efron and Morris 1975). S_theta is estimated as the excess of the weighted scatter of y over v (0 if negative). The predicted fall of a place is B (y - ybar): the size of the regression to the mean a hot-spot programme should expect without doing anything. That B* lies in [0, 1], minimises the loss and never does worse than no shrinkage is proved in Lean: see proof("ebshr").

Usage

ebshr(y, v = NULL, w = NULL, b = NULL)

Arguments

y

The first-period counts (or rates): finite numbers, at least two.

v

NULL for the weighted mean of y (Poisson noise), or the noise variance, a number from 0.

w

NULL, or one weight from 0 per place, positive sum.

b

NULL for B*, or a shrinkage factor from 0 to 1.

Value

A list of class ebshr: fit, a named vector: mean, total (the weighted scatter of y), noise (v), signal (the estimated variance of the rates), bstar and b (the factor used); est, a data frame with one row per place: y, shrunk and fall.

References

Efron, B. and Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association, 70(350), 311-319. doi:10.1080/01621459.1975.10479864

Examples

y <- c(40, 12, 9, 25, 7, 31, 5, 18)
s <- ebshr(y)
s
#> ebshr: shrinkage 0.1327 (B* 0.1327) toward 18.38; signal variance 120.1, noise 18.38
#>   y    shrunk        fall
#>  40 37.130656  2.86934447
#>  12 12.845876 -0.84587611
#>   9 10.243935 -1.24393546
#>  25 24.120952  0.87904773
#>   7  8.509308 -1.50930836
#>  31 29.324834  1.67516642
#>   5  6.774681 -1.77468126
#>  18 18.049757 -0.04975742
s$est[which.max(y), ]
#>    y   shrunk     fall
#> 1 40 37.13066 2.869344