---
title: "Small-area rates: expected counts, shrinkage, funnels and spatial autocorrelation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Small-area rates: expected counts, shrinkage, funnels and spatial autocorrelation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# the vignette runs in full in GitHub CI and for the package tarball (NOT_CRAN), on r-universe
# (MY_UNIVERSE), on the package website (IN_PKGDOWN) and in an interactive session; when CRAN
# re-builds it in its check, none of these is set and only this chunk runs
if (!interactive() && !any(nzchar(Sys.getenv(c("NOT_CRAN", "MY_UNIVERSE", "IN_PKGDOWN"))))) {
  knitr::knit_exit()
}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.8)
library(rmoriebricklayer)
```

Published tables often give a count per region, division or facility. Comparing those
counts fairly needs three steps: adjust for the population each unit serves (and its
composition), stabilise the ratios of small units, and judge each unit against limits
that widen as its expected count shrinks. This vignette walks the four tools that do that.

## Expected counts by indirect standardisation

If the overall rate applied everywhere, how many events would each unit have? That is the
*expected count*. Adjusting for size alone:

```{r small-area-rates-expected}
d <- data.frame(
  region = rep(c("North", "South", "East"), each = 2),
  age = rep(c("18 to 24", "25 to 49"), 3),
  n = c(30, 45, 12, 60, 8, 20),
  pop = c(1000, 6000, 900, 9000, 400, 3500)
)
expcnts(d$n, d$pop, d$region)
```

Adjusting for size **and** age composition is indirect standardisation: each stratum's
overall rate is applied to the unit's population in that stratum, then summed. This is the
comparison that says something about the unit rather than about who lives there:

```{r small-area-rates-expected-strata}
expcnts(d$n, d$pop, d$region, strata = d$age)
```

## Direct standardisation, for populations of different ages

Indirect standardisation asks what each unit would expect at the overall rates. Direct
standardisation asks what rate each population would have with one standard age structure. It
is the comparison to use when the populations' ages differ a great deal. First Nations, Metis and
Inuit populations, for example, are younger than the non-Indigenous population.
`agestd()` gives the rates with the Fay-Feuer gamma interval, which holds for small
counts. `rategap()` compares two of them. The counts below are made up, to show the computation:

```{r small-area-rates-direct}
r <- agestd(events = c(4, 9, 11, 6, 3, 40, 70, 90, 80, 60),
                     pop = c(9000, 8000, 6000, 4000, 2000, 300000, 320000, 350000, 340000, 300000),
                     std = c(30, 25, 20, 15, 10),
                     grp = rep(c("Inuit", "Non-Indigenous"), each = 5))
r
rategap(r, "Inuit", "Non-Indigenous")
```

For survey percentages rather than counts, `stdprev()` does the same.

## Standardised ratios with exact intervals

The ratio of observed to expected is the standardised incidence ratio (SIR). `stdinc()` gives
the exact Poisson interval, which is what makes a small count honest:

```{r small-area-rates-sir}
stdinc(obs = c(30, 12, 3), expd = c(20, 14, 5),
    area = c("North", "South", "East"))
```

An observed count of three carries almost no information, and the interval says so rather
than hiding it:

```{r small-area-rates-sir3}
stdinc(3, 5)
```

## Empirical Bayes shrinkage

A tiny unit's raw ratio is extreme by chance. Empirical Bayes shrinks each ratio toward the
overall experience in proportion to how little its own data say: the pooled ratios share a
Gamma prior whose parameters are estimated from all units, and each unit's posterior mean
is a weighted average of its own ratio and the prior mean, weighted by its expected count
(Clayton and Kaldor 1987).

```{r small-area-rates-eb}
ebrates(obs = c(30, 45, 2), expd = c(25, 50, 0.5),
         area = c("North", "South", "Tiny"))
```

The tiny area's raw ratio of 4 is pulled sharply toward 1; the large areas barely move.

## Funnel plots

A funnel plot draws control limits that narrow as the expected count grows (Spiegelhalter
2005): a ratio of 2 means nothing at an expected count of 2 and a great deal at 200. Units
outside the limits are the ones worth a second look; the rest are consistent with the
common rate.

```{r small-area-rates-funnel}
fnllim(c(2, 20, 200))
```

```{r small-area-rates-funnel-plot}
e <- seq(2, 200, length.out = 60)
lim <- fnllim(e)
plot(e, lim$upper_95, type = "l", ylim = c(0, 3), xlab = "expected count",
     ylab = "observed / expected", main = "Funnel limits")
lines(e, lim$lower_95)
lines(e, lim$upper_998, lty = 2)
lines(e, lim$lower_998, lty = 2)
abline(h = 1, col = "grey")
points(c(20, 14, 5), c(30, 12, 3) / c(20, 14, 5), pch = 19, col = 2)
```

The `levels` argument sets which limits are drawn; 95% and 99.8% (about two and three
standard deviations) are the conventional pair.

## Spatial autocorrelation

Rates of neighbouring areas are rarely independent. Moran's *I* measures whether similar
values cluster in space, against a permutation null; `neighbours` is a list giving each
area's adjacent areas:

```{r small-area-rates-moran}
nb <- list(2L, c(1L, 3L), c(2L, 4L), c(3L, 5L), c(4L, 6L), 5L)   # six areas in a line

set.seed(1)
morani(c(1, 2, 3, 4, 5, 6), nb, nperm = 999L)[c("I", "p_value")]   # a gradient

set.seed(1)
morani(c(1, 6, 1, 6, 1, 6), nb, nperm = 999L)[c("I", "p_value")]   # alternating
```

The null expectation of *I* is not zero but $-1/(n-1)$:

```{r small-area-rates-moran-null}
-1 / (6 - 1)
```

Positive autocorrelation in the residual ratios means neighbouring units share something
the standardisation did not capture, and a funnel that treats them as independent is
optimistic.

## From a published table to a judgement

1. `expcnts()` with the strata the table publishes.
2. `stdinc()` for the ratios and their exact intervals.
3. `ebrates()` where small units dominate.
4. `fnllim()` to see which units are outside the common experience.
5. `morani()` on the residuals before trusting independence.

## References

Clayton, Kaldor (1987). Empirical Bayes estimates of age-standardized relative risks for
use in disease mapping. *Biometrics* 43(3).
Spiegelhalter (2005). Funnel plots for comparing institutional performance. *Statistics in
Medicine* 24(8).
Moran (1950). Notes on continuous stochastic phenomena. *Biometrika* 37(1/2).