# Evaluation of four gamma-based methods for calculating confidence intervals for age-adjusted mortality rates when data are sparse – Population Health Metrics

May 7, 2022

### Technical appendix

The number of deaths D reported for a given area and time period is assumed to be Poisson-distributed, with mean ({mathbb{E}})(D) and variance ({mathbb{V}})(D) satisfying the equality ({mathbb{E}})(D) = λP = ({mathbb{V}})(D), where P denotes the population denominator [1]. The age-specific or crude death rate R, defined as the ratio D/P, is usually multiplied by 100,000 and reported as a rate per 100,000 population.

### Poisson-gamma relationship

For a positive integer x ≤ P, it can be shown [2] that there exists a gamma random variable G such that﻿ ({mathbb{E}})(G) = x =﻿ ({mathbb{V}})(G) and

$$Pr left( {D ge x|lambda } right) = Pr left( {G le lambda P|x} right)$$

(A1)

Recall that if G is gamma-distributed with shape parameter α > 0 and scale parameter β > 0, then its mean and variance are ({mathbb{E}})(G) = αβ and ({mathbb{V}})(G) = αβ2. Conversely, the parameters are given by α = ﻿({mathbb{E}})(G)2/({mathbb{V}})(G) and β = ({mathbb{V}})﻿(G)/({mathbb{E}})(G). Thus, with ﻿({mathbb{E}})(G) = x = ({mathbb{V}})(G) in Eq. A1, the corresponding gamma distribution has α = x and β = 1.

For the rate R = D/P, with P = p, y = x/p, and v = x/p2, Eq. A1 becomes

$$Pr (R ge y|lambda ) = Pr (Z le lambda |y,v)$$

(A2)

where Z = G/P is gamma-distributed with mean y and variance v.

### Gamma CI for age-specific and crude rates

When D = x is observed, the ratio y = x/p is an estimate of ({mathbb{E}})(R) = λ. An equal-tailed 100(1 − a) percent CI [L(y), U(y)] for the parameter λ, e.g., with a = 0.05, is obtained as a solution to the following two equations:

$$Pr left[ {R ge y|lambda = Lleft( y right)} right] = a/2$$

(A3a)

$$Pr left[ {R le y|lambda = Uleft( y right)} right] = a/2$$

(A3b)

Eqs. A3a and A3b follow from looking upon L(y) as the largest λ for which Pr(R ≥ y|λ) ≤ a/2 and U(y) as the smallest λ for which Pr(R ≤ y|λ) ≤ a/2; see [22] and theorem 9.2.3.a in [23].

From Eqs. A2 and A3a,

$$a/2 = Pr [R ge y|lambda = L(y)] = Pr [Z le Lleft( y right)|y,v]$$

where Z is gamma-distributed with mean y and variance v, i.e., with parameters α = y2/v = x and β = v/y = 1/p. Thus, the lower CI limit L(y) is obtained as the (a/2)-quantile of the gamma(x, 1/p) distribution. For y = 0 = x, L(0) = 0 by convention.

Similarly, from Eqs. A2 and A3b,

$$1 – left( {a/2} right) = Pr [R > y|lambda = U(y)] quadquadquadquad, = Pr [R ge y + 1/p|lambda = U(y)] quadquadquadquad, = Pr [Z^{prime} le U(y)|y^{prime},v^{prime}]$$

where the second equality is due to x being a positive integer, so that D/p > x/p if and only if D/p ≥ (x + 1)/p, and Z′ is a gamma random variable with mean y′ = y + 1/p and variance v′ = v + 1/p2. Because y′2/v′ = x + 1 and v′/y′ = 1/p, the upper CI limit U(y) is obtained as the (1 − a/2)-quantile of the gamma(x + 1, 1/p) distribution.

### Approximate gamma CIs for age-adjusted rates

With n age groups, let Di denote the number of deaths for group i. The Di are assumed to be independent Poisson random variables, and the age-specific rates Ri are defined as the ratios Di/Pi, with means ({mathbb{E}})(Ri) = λi and variances ({mathbb{V}})(Ri) = λi/pi.

Let πi denote the size of group i in the reference population, e.g., the projected year 2000 US population [24]. Let wi denote the relative proportions for group i in the reference population: wi = πi/∑πj. The age-adjusted death rate R′ is defined as

$$R^{prime} = sum w_{i} R_{i} = sum left( {w_{i} /P_{i} } right)D_{i}$$

Given the parameters λi and denominators Pi = pi, the age-adjusted rate R′ has mean ﻿({mathbb{E}})(R′) = λ′ = ∑wi λi and variance ﻿({mathbb{V}})(R′) = ∑wi2 λi/pi.

Fay–Feuer interval. Fay and Feuer [7] assume that Eq. A2 holds approximately for the age-adjusted rate R′, so that, for y = ∑(wi/pi) xi and v = ∑(wi/pi)2 xi,

$$Pr (R^{prime} ge y|lambda^{prime}) approx Pr (Z le lambda^{prime}|y,v)$$

(A4)

where Z is gamma-distributed with ({mathbb{E}})(Z) = y and ({mathbb{V}})﻿(Z) = v, i.e., with α = y2/v and β = v/y. As for the crude rate R, an equal-tailed 100(1 − a) percent CI for λ′ solves the equations:

$$Pr left[ {R^{prime} ge y|lambda^{prime} = Lleft( y right)} right] = a/2$$

(A5a)

$$Pr left[ {R^{prime} le y|lambda^{prime} = Uleft( y right)} right] = a/2$$

(A5b)

From Eqs. A4 and A5a, the lower limit L(y) can be resolved approximately from the lower tail probability of a gamma distribution with parameters α = y2/v and β = v/y, again with the convention that L(0) = 0.

For the upper bound, note that a unit increment in the observed number of deaths xj within group j results in the addition of the quantity wj/pj to the age-adjusted rate y = ∑(wi/pi) xi. Because such a unit increment could be realized in any of the n groups,

$$Pr [R^{prime} > y|lambda^{prime} = U(y)] ge Pr [R^{prime} ge y + kappa_{0} |lambda^{prime} = U(y)]$$

where κ0 = max{wj/pj}. From Eq. A4, the right-hand side in this last inequality is approximately equal to Pr[Z′ ≤ U(y)|y′, v′], where Z′ is gamma-distributed with mean y′ = y + κ0 and variance v′ = v + κ02. Thus, an upper CI limit U(y) can be resolved from the upper tail probability of a gamma distribution with shape parameter α = y′2/v′ and scale parameter β = v′/y′. Fay and Feuer [7] make the conjecture that the approximate gamma CI thus constructed remains conservative. Although this conjecture remains unproven, findings from the many simulation studies to date continue to support it, e.g., [4,5,6].

Tiwari modification. Tiwari et al. [8] developed a modification to the Fay–Feuer method described above by distributing an average increment 1/n uniformly across all age groups instead of a unit increment in a single age group:

$$y^{prime} = sum left( {w_{i} /p_{i} } right)left( {x_{i} + frac{1}{n}} right) = y + frac{1}{n}sum w_{i} /p_{i}$$

Thus, with κ1 = n−1 ∑wi/pi and κ2 = n−1 ∑(wi/pi)2, the gamma random variable Z′ above now has mean y′ = y + κ1 and variance v′ = v + κ2. The Tiwari modification reduces the CI width relative to the Fay–Feuer method, because

$$Pr [R^{prime} ge y + kappa_{1} |lambda^{prime} = U(y)] ge Pr [R^{prime} ge y + kappa_{0} |lambda^{prime} = U(y)]$$

However, the resulting CI sometimes fails to retain the nominal coverage level, e.g., [4].

Fay–Kim modification. Fay and Kim [10] more recently developed a mid-p version of the Fay–Feuer CI. A modification of exact CIs from discrete data, mid-p CIs trade-off guaranteed nominal coverage in all of the parameter space (which tends to result in overly wide CIs) for proximity to nominal coverage (and narrower CIs) for most parameter values.

For the mid-p interval, a solution to the following equations is sought:

$$Pr left[ {R^{prime} > y|lambda^{prime} = Lleft( y right)} right] quadquadquadquad + left( {1/2} right) times Pr left[ {R^{prime} = y|lambda^{prime} = Lleft( y right)} right] = a/2$$

(A7a)

$$Pr left[ {R^{prime} < y|lambda^{prime} = Uleft( y right)} right] quadquadquadquad + left( {1/2} right) times Pr left[ {R^{prime} = y|lambda^{prime} = Uleft( y right)} right] = a/2$$

(A7b)

Drawing B = b from a Bernoulli distribution with Pr(B = 1) = 1/2, Fay and Kim [10] define the mid-p version of the Fay–Feuer CI using the following gamma distribution:

$${text{gamma}}_{{text{mid-p}}} = b times {text{gamma}}left( {y^{2} /v,v/y} right) quadquadquadquadquadquad + left( {1 – b} right) times {text{gamma}}left( {y^{{prime}{2}} /v^{prime},v^{prime}/y^{prime}} right)$$

where y′ = y + κ0 and v′ = v + κ02 are as in the Fay–Feuer construction, above. Thus, the lower and upper limits are defined as the (a/2)th and (1−a/2)th quantiles of gamma mid-p. The special case y = 0 is addressed using L(0) = 0 and U(0) defined as the (1−a)th quantile of the gamma(y2/v, v/y) distribution. R syntax is provided to solve for L(y) and U(y) numerically [10].

$$Pr left( {D_{{{text{adj}}}} ge x_{{{text{adj}}}} |lambda^{prime}} right) = Pr left( {G_{{{text{adj}}}} le lambda^{prime}p_{{{text{adj}}}} |x_{{{text{adj}}}} } right)$$

(A6)

Because y2/v will generally not be integer, xadj is defined as the nearest integer instead (although this is not strictly necessary), and the equality in Eq. A6 is assumed to hold approximately. Either way, one proceeds as for the crude rate to derive CI limits L(y) and U(y) for λ′ as the (a/2)-quantile of the gamma(xadj, 1/padj) distribution and the (1 − a/2)-quantile of the gamma(xadj + 1, 1/padj), respectively.

Exact intervals. When there is a constant scalar c > 0 such that pi = i for all i, the age-adjusted rate equals the overall crude rate, and the above CIs reduce to the exact gamma CI for λ = p1λi pi where p = ∑pi and the total number of deaths D = ∑Di follows a Poisson distribution with mean λp. In particular, when y = 0, v = 0 and xadj is undefined. However, because the age-adjusted rate equals the crude rate in this case, the limits of all three approximate gamma CIs for the age-adjusted rate are defined to be those of the exact gamma CI for the crude rate, with p = ∑pi and x = ∑xi = 0. Thus, in this extreme case, L(0) = 0 and U(0) is the (1 − a/2)-quantile of the gamma(1, 1/p) distribution.

Anderson–Rosenberg CI as a modification of the Fay–Feuer CI. The Anderson–Rosenberg construction can be seen to follow that of the Fay–Feuer CI, with a gamma-distributed Z′′ that has mean y′′ = y + κ and variance v′′ = v + κ2, where κ = κ3 = 1/padj instead of κ = κ0 = max{wj/pj}. Indeed, with 1/padj = v/y and xadj = y2/v,

$$frac{{v^{primeprime} }}{{y^{primeprime} }} = frac{{v(y^{2} + v)/y^{2} }}{{(y^{2} + v)/y}} = frac{v}{y} = frac{1}{{p_{{{text{adj}}}} }}$$

$${text{and}}quad frac{{y^{{primeprime}{2}} }}{{v^{primeprime}}} = frac{{(y^{2} + v)^{2} /y^{2} }}{{v(y^{2} + v)/y^{2} }} = frac{{y^{2} + v}}{v} = x_{{{text{adj}}}} + 1$$

Furthermore, 1/padj can be expressed as follows:

begin{aligned}frac{1}{{p_{{{text{adj}}}} }} &= frac{v}{y} = sum left( {w_{i} /p_{i} } right)xi_{i} quad{text{with}}\xi_{i} &= frac{{(w_{i} /p_{i} )x_{i} }}{{sum (w_{j} /p_{j} )x_{j} }};{text{and}};sum xi_{i} = 1.end{aligned}

As a result,

$$y^{primeprime} = y + frac{1}{{p_{{{text{adj}}}} }} = sum (w_{i} /p_{i} )(x_{i} + xi_{i} )$$

and ﻿the Anderson–Rosenberg method is seen to result in incrementing the age-specific death counts from xi to xi + ξi, whereas in the Fay–Feuer method only the count xi* for the age group i* for which wi*/pi* = max{wj/pj} is incremented—and in the Tiwari modification, the age-specific counts are incremented from xi to xi + ζi, where ζi = 1/n. Additionally,

$$kappa_{3} = frac{1}{{p_{{{text{adj}}}} }} = sum (w_{i} /p_{i} )xi_{i} le max { w_{i} /p_{i} } = kappa_{0}$$

since ﻿∑ξi = 1. Thus, like the Tiwari modification, the Anderson–Rosenberg construction reduces the CI width relative to the Fay–Feuer method:

$$Pr [R^{prime} ge y + kappa_{3} |lambda^{prime} = U(y)] ge Pr [R^{prime} ge y + kappa_{0} |lambda^{prime} = U(y)]$$

Two questions emerge from the above derivations:

1. (1)

Under what circumstances does the Anderson–Rosenberg method result in a shorter CI than the Fay–Feuer method that retains nominal coverage?

2. (2)

Since both the Anderson–Rosenberg and Tiwari methods result in narrower CIs than the Fay–Feuer method, when is one preferable to the other?

To partially answer question 2, note that the Anderson–Rosenberg CI would be narrower than the Tiwari CI if (but not only if) κ3 ≤ κ1, as that ensures

$$Pr [R^{prime} ge y + kappa_{3} |lambda^{prime} = U(y)] ge Pr [R^{prime} ge y + kappa_{1} |lambda^{prime} = U(y)]$$

By definition, the condition κ3 ≤ κ1 is realized when

$$sum (w_{i} /p_{i} )xi_{i} le frac{1}{n}sum (w_{i} /p_{i} )$$

﻿which is equivalent to

$$frac{1}{n}sum (w_{i} /p_{i} )^{2} x_{i} le left{ {frac{1}{n}sum (w_{i} /p_{i} )x_{i} } right} left{ {frac{1}{n}sum (w_{i} /p_{i} )} right}$$

This last condition indicates that the slope of the line from the simple regression of the weight-adjusted age-specific death rates wi (xi/pi) = (wi/pi) xi onto the weights wi/pi is negative or zero. This could be verified upfront for any set of age-adjustment weights (w1, …, wn) and population distribution (p1, …, pn), and it would be sufficient to ensure that the Anderson–Rosenberg CI will be narrower than the Tiwari CI. Of course, this leaves the issue of efficiency unresolved, as it would theoretically be possible for either the Anderson–Rosenberg or the Tiwari CIs to be so narrow as to fail to retain nominal coverage. The empirical simulations investigate situations where this may occur. In addition, these two CI methods are compared to the more recent Fay–Kim mid-p modification.