An easy approach to the Robins-Breslow-Greenland variance estimator

Abstract

The Mantel-Haenszel estimate for the odds ratio (and its logarithm) in stratified case control studies lacked a generally acceptable variance estimate for many years. The Robins-Breslow-Greenland estimate has met this need, but standard textbooks still do not provide an explanation of how it is derived. This article provides an accessible derivation which demonstrates the link between the Robins-Breslow-Greenland estimate and the familiar Woolf estimate for the variance of the log odds ratio, and which could easily be included in Masters level courses in epidemiology. The relationships to the unconditional and conditional maximum likelihood estimates are also reviewed.

Introduction

The Mantel-Haenszel (MH) estimate for the summary odds ratio across several 2 × 2 tables, ψ _MH , was proposed in 1959 [1]. Over twenty years later the lack of a robust estimate for its variance was still being noted [2], yet only a few years afterwards, Robins, Breslow and Greenland introduced their now generally-accepted variance estimator [3] for the Mantel-Haenszel log-odds ratio (denoted by the RBG estimate). This replaced estimation of confidence limits based on the unsatisfactory test-based procedure of Miettinen or the computationally intensive Cornfield type limits which had hitherto been used.

While a useful review of Mantel-Haenszel methods has been published, including some aspects of the historical development towards the RBG estimator [4] the formal derivations by Robins, Breslow & Greenland [3] and Phillips & Holland [5] are not, in the view of this author, easily comprehended. The former omits steps in the argument, while the latter appeals to descending factorial powers. Possibly it is no surprise that even modern textbooks [6, 7] merely state the RBG formula without deriving it.

While other variance estimators exist, some are ad hoc, such as the application of the cohort study formula to case-control data suggested by Clayton and Hills [8], only apply to the large few strata case [9] or are closely related to the RBG estimator [10]. One rather different exception is Sato's formula [11] but this procedure gives confidence limits directly in the odds ratio scale.

It is the intention of this article to present an informal derivation of the RBG estimator as an extension of the familiar variance formula of Woolf [12], and which could readily be included in standard textbooks of epidemiology or biostatistics. I will describe this from the perspective of a case-control study.

Analysis

How does the Mantel-Haenszel estimate arise?

Consider a stratified case-control study for which the i^th of k independent tables is:

Neglecting constants, the unconditional likelihood for the i^th table is:

where in the i^th table θ _i = probability of exposure if a case and ø _i = probability of exposure if a control.

The maximum likelihood estimate (MLE) for ø _i is given by b _i /(b _i + d _i ) and if we re-parameterise θ _i as ψ ø _i / [ψ ø _i + (1 - ø _i )], where ψ is the odds ratio (assumed common to all tables), the contribution to the overall log likelihood made by terms involving ψ is:

∑ a _i ln {ψ ø _i /[ψ ø _i + (1 - ø _i )]} + c _i ln {1/[ψ ø _i + (1 - ø _i )]}.

Differentiating with respect to ψ and equating to zero, and rearranging (noting that a _i + c _i = n _1i )

we obtain:

∑ a _i ln {ψ ø _i /[ψ ø _i + (1 - ø _i )]} + c _i ln {1/[ψ ø _i + (1 - ø _i )]}.

i.e., ∑ {a _i - n _1i ψ ø _i /[ψ ø _i + (1 - ø _i )]} = 0

i.e., ∑ {[ψ a _i ø _i + a _i - a _i ø _i - n _1i ψ ø _i ]/[ψ ø _i + (1 - ø _i )]} = 0.

This must be solved numerically to obtain the MLE for ψ, but if the denominators do not vary too much across the tables we merely have to solve:

∑ [ψ a _i ø _i + a _i - a _i ø _i - n _1i ψ ø _i ] = 0

i.e., ∑ [ψ(a _i - n _1i )ø _i + a _i (1 - ø _i )] = 0

or, ∑ a _i (1 - ø _i ) = ∑ ψ(n _1i - a _i )ø _i

giving, ∑ a _i (1 - ø _i ) = ψ ∑ (n _1i - a _i )ø _i

and since, ψ _i = b _i /(b _i + d _i ) = b _i /n _0i

This can be used as a first approximation to find the MLE (if there is only one table then ψ is the unconditional MLE = ad/bc). Now in stratified case-control studies with a constant ratio, r, of controls to cases, the total number of subjects in each stratum is given by n _i = n _0i (1 + r), so n _0i = n _i /(1 + r). A constant r will be achieved by design if there is caliper matching; otherwise – as with a post-stratified analysis – this will be only approximately true. The term (1 + r) can then be cancelled and we are left with:

The MH estimator is therefore a first approximation to the unconditional MLE in the large strata case with a constant control:case ratio across strata. However the MH estimator actually coincides with the conditional MLE for the matched pairs design, as outlined, for example on page 164 of Breslow & Day [2].

The sensitivity to variation in the ø _i and constancy of the control:case ratio is not high, as shown by the data in Table 1. In a sense this would be expected because for the most sparse (e.g., pair-matched) data the control:case ratio will be constant, and while the ø _i then have maximum variance – being only 0 and 1, the MH estimate coincides with the conditional MLE. Conversely, for large strata the control:case ratio will vary, but the variance of the ø _i will be less and the MH estimate will then approximate the unconditional MLE.

Table 1 Simulated case-control data with true odds ratio = 5

Deriving the variance of the Mantel-Haenszel estimate

Consider again the i^th 2 × 2 table, giving the frequencies in each cell:

For odds ratio , estimated for a single table by the cross-product ratio a _i d _i /b _i c _i , application of the delta method gives Woolf's logit-based formula [8]:

with n _i = a _i + b _i + c _i + d _i and,

The delta method is a widely used procedure in statistics when an approximation is needed for the variance of a function of a variable whose variance is known. In this instance the variable with known variance is a proportion p, and the function is the logit. The basic delta method formula is: var(y) ≈ (dy/dx)² var(x) from which, if y = logit(p = ln[p/(1 - p)],

var(y) ≈ (1/p + 1/(1 - p))² p(1 - p)/n

= (1/p + 1/(1 - p))1/n

= (1/a + 1/b)

if p = a/n and n = a + b.

Here we have two independent proportions (the proportion of cases and controls exposed) and Woolf's formula is obtained by estimating the variances of the separate logits and adding them.

For k such 2 × 2 tables, each representing a separate stratum, the Mantel-Haenszel pooled estimate of the common odds ratio ψ is given by:

Hence ψ _MH is a weighted average of the stratum-specific odds ratios. The weights approximate the inverse of the variance of each _i if the true value of ψ = 1. Note that the assumption here of a common odds ratio is not required for the Mantel-Haenszel test.

To derive the variance, in addition to the approximation involved in application of the delta rule, an assumption is also made that each stratum-specific odds ratio is close enough to the Mantel-Haenszel pooled estimate to permit terms like a _i d _i /b _i c _i to be replaced by ψ _MH .

We then proceed by obtaining an approximation which avoids zeros in the formula for var[ln()]. The motivation for this can be seen by comparing the weights for ψ _MH – which are unaffected by zeros except for deleting such strata – whereas if Woolf's variances were used, the result would be indeterminate if cells with zeros were present.

Taking the weights as constant,

Assuming a common odds ratio ψ, estimated by ψ _MH , this can be written as:

Leading to a formula suggested by Hauck [9]:

As mentioned above, a problem with this formula is that it fails if cell entries are zero. However we can proceed further by re-writing the formula as:

On substituting 1/ψ _MH for (b _i c _i /a _i d _i ):

Now if the rows of the 2 × 2 table are interchanged, the variance stays the same. But a similar argument to that above leads to:

(Note that the new odds ratio formed by exchanging rows is just 1/ψ _MH .) "The" variance, V, of ln(ψ _MH ) is therefore taken to be the mean of the two estimates [13] as follows:

Let R = ∑ (a _i d _i /n _i ) and S = ∑ (b _i c _i /n _i ). On substituting into the two variance formulae:

Next, divide the top and bottom by S² and move the term outside the brackets to obtain:

which is eq. 9 in Phillips & Holland [5].

If we now put

P_i = (a _i + d _i )/n _i and Q_i = (b _i + c _i )/n _i with R_i = a _i d _i /n _i and S_i = b _i c _i /n _i

then

which on multiplying out the brackets, rearranging and noting that R/S = ψ _MH , gives:

This is the RBG formula!

When there is only one stratum, this reduces to (1/a + 1/b + 1/c + 1/d) which is the familiar logit based formula of Woolf and which approaches 0 as the sample size increases, assuming a finite true odds ratio. Clearly as the RBG variance estimate is a finite sum of such estimators the RBG estimate will also approach 0, for large strata.

The RBG estimator was derived above on the assumption that the stratum-specific odds ratio estimates could be liberally replaced by the common value, in turn estimated by ψ _MH ; both assumptions are reasonable with large samples per stratum. However, the success of the RBG formula derives from its being applicable also to the sparse data case.

To see this, consider a matched-pairs case control study. The capital letters denote the frequency of case-control pairs.

In such a study each stratum has only two observations. The table can be decomposed into four types of "unmatched" table according to the exposure category of the case and the control, the frequency of each type being given by the frequency of the corresponding case-control pairs:

Only the B such tables with a _i = d _i = 1 and the C such tables with b _i = c _i = 1 contribute to the estimate of the odds ratio. Note that these are disjoint sets of tables.

Under these circumstances: ψ _MH = B/C which coincides with the conditional MLE and:

a) the middle term of the RBG formula vanishes because if b _i c _i = 1 then (a _i + d _i ) = 0, and if a _i d _i = 1 then (b _i + c _i ) = 0

b) R = ∑ a _i d _i /n _i = B/2 & S = ∑ b _i c _i /n _i = C/2

c) There are B terms in which a _i d _i (a _i + d _i ) = 2 C terms in which b _i c _i (b _i + c _i ) = 2

giving:

V = B/B² + C/C² = 1/B + 1/C

This is not only the familiar logit based formula for the variance of the log odds ratio for matched pairs, but is also the variance of the conditional maximum likelihood estimate. This is asymptotically consistent from the general properties of a MLE (and it's easy to see that as the number of tables increases, V → 0).

In other words, the RBG formula, though derived here without assuming validity in the sparse case, does in fact possess this property.

Table 1 shows how closely the conditional maximum likelihood estimate, unconditional maximum likelihood estimate, and MH estimate agree, despite varying ø _i and control:case ratio.

Conclusion

The Mantel-Haenszel estimate of the odds ratio approximates the maximum likelihood estimate for large, few strata and coincides with the conditional maximum likelihood estimate for the sparse data (matched pairs) case.

The RBG formula is the estimator of choice for the variance of the Mantel-Haenszel log-odds-ratio because it applies both in the large few strata case and in the many sparse strata case (as in matched pairs analysis), when the RBG variance estimate actually coincides with the conditional maximum likelihood variance estimate.

Moreover the RBG formula reduces to familiar standard forms for a single stratum and for matched pairs.

Formal derivation of the RBG formula is tricky but an informal, accessible derivation is possible as outlined above, which uses nothing more advanced than the delta method for approximating a variance.