Working through a confounded analysis

The data

In 1972–74, a survey in Whickham, England recorded whether each participant smoked. Twenty years later the survivors were identified. The data below are the 1314 women in that cohort, and they are a standard illustration of Simpson’s paradox (Appleton, French and Vanderpump, 1996, The American Statistician 50, 340–341).

Start where a student would, ignoring everything except smoking and death.

crude <- epi2x2(139, 443, 230, 502,
                exposure = c("Smoker", "Non-smoker"),
                outcome  = c("Dead", "Alive"))
crude
#>             Dead  Alive  Total
#> Smoker       139    443    582
#> Non-smoker   230    502    732
#> Total        369    945   1314

The crude analysis

odds_ratio(crude)

Odds ratio

Step 1. Odds of Dead among the exposed (Smoker)

\[ \begin{aligned} \text{odds}_1 &= \frac{a}{b} \\ &= 139 / 443 \\ &= 0.3138 \end{aligned} \]

Step 2. Odds of Dead among the unexposed (Non-smoker)

\[ \begin{aligned} \text{odds}_0 &= \frac{c}{d} \\ &= 230 / 502 \\ &= 0.4582 \end{aligned} \]

Step 3. Odds ratio

\[ \begin{aligned} OR &= \frac{\text{odds}_1}{\text{odds}_0} \\ &= 0.3138 / 0.4582 \\ &= 0.6848 \end{aligned} \]

Equivalently OR = ad / bc, which is why the odds ratio is the same whether you condition on exposure or on outcome. That symmetry is what makes it usable in case-control studies.

Step 4. Standard error of log(OR), Woolf’s method

\[ \begin{aligned} SE &= \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \\ &= \operatorname{sqrt}(1/139 + 1/443 + 1/230 + 1/502) \\ &= 0.1257 \end{aligned} \]

Step 5. 95% confidence interval, lower limit

\[ \begin{aligned} \text{lower} &= \exp\left(\log OR - z \cdot SE\right) \\ &= \exp(-0.3786 - 1.96 \times 0.1257) \\ &= 0.5353 \end{aligned} \]

Step 6. 95% confidence interval, upper limit

\[ \begin{aligned} \text{upper} &= \exp\left(\log OR + z \cdot SE\right) \\ &= \exp(-0.3786 + 1.96 \times 0.1257) \\ &= 0.8761 \end{aligned} \]

Result: \(OR = 0.6848\) (95% CI 0.5353 to 0.8761)

• These data are tabulated as a cohort, where risk is estimable. The risk ratio here is 0.7601 against an odds ratio of 0.6848 – the odds ratio is the more extreme of the two, and always will be. The two converge only when the outcome is rare; baseline risk here is 0.3142.

Three things are worth noticing about that output.

The odds ratio is 0.68, and the confidence interval excludes 1. Read naively, smoking is protective, and significantly so.

The derivation shows the odds in each group separately before dividing them. A student who computes 139/443 = 0.31 and stops has done step 1 correctly, and telling them “wrong, the answer is 0.68” hides that from both of you.

The note at the bottom appeared without being asked for. Because the data are tabulated as a cohort, risk is estimable, so the function reports what the risk ratio would be and reminds you that the odds ratio is always the more extreme of the two.

Since risks are estimable here, compute them directly. There is no need to see every step again, so turn the detail down:

options(epibyhand.verbose = 1)
risk_ratio(crude)

Risk ratio

Step 1. Risk of Dead among the exposed (Smoker)

\[ \begin{aligned} R_1 &= \frac{a}{a + b} \\ &= 0.2388 \end{aligned} \]

Step 2. Risk of Dead among the unexposed (Non-smoker)

\[ \begin{aligned} R_0 &= \frac{c}{c + d} \\ &= 0.3142 \end{aligned} \]

Step 3. Risk ratio

\[ \begin{aligned} RR &= \frac{R_1}{R_0} \\ &= 0.7601 \end{aligned} \]

Step 4. Standard error of log(RR)

\[ \begin{aligned} SE &= \sqrt{\frac{1}{a} - \frac{1}{a+b} + \frac{1}{c} - \frac{1}{c+d}} \\ &= 0.092 \end{aligned} \]

The interval is built on the log scale because RR is a ratio: its sampling distribution is skewed, but log(RR) is roughly normal.

Step 5. 95% confidence interval, lower limit

\[ \begin{aligned} \text{lower} &= \exp\left(\log RR - z \cdot SE\right) \\ &= 0.6347 \end{aligned} \]

Step 6. 95% confidence interval, upper limit

\[ \begin{aligned} \text{upper} &= \exp\left(\log RR + z \cdot SE\right) \\ &= 0.9102 \end{aligned} \]

Result: \(RR = 0.7601\) (95% CI 0.6347 to 0.9102)

• A risk ratio requires that everyone was followed for the same period. If follow-up time varied, use a rate ratio on person-time instead.

risk_difference(crude)

Risk difference

Step 1. Risk of Dead among the exposed (Smoker)

\[ \begin{aligned} R_1 &= \frac{a}{a + b} \\ &= 0.2388 \end{aligned} \]

Step 2. Risk of Dead among the unexposed (Non-smoker)

\[ \begin{aligned} R_0 &= \frac{c}{c + d} \\ &= 0.3142 \end{aligned} \]

Step 3. Risk difference

\[ \begin{aligned} RD &= R_1 - R_0 \\ &= -0.0754 \end{aligned} \]

Step 4. Standard error of RD

\[ \begin{aligned} SE &= \sqrt{\frac{R_1(1-R_1)}{a+b} + \frac{R_0(1-R_0)}{c+d}} \\ &= 0.0246 \end{aligned} \]

No log transform here. A difference can be negative, so it is already on a scale where the normal approximation applies directly.

Step 5. 95% confidence interval, lower limit

\[ \begin{aligned} \text{lower} &= RD - z \cdot SE \\ &= -0.1237 \end{aligned} \]

Step 6. 95% confidence interval, upper limit

\[ \begin{aligned} \text{upper} &= RD + z \cdot SE \\ &= -0.0271 \end{aligned} \]

Result: \(RD = -0.0754\) (95% CI -0.1237 to -0.0271)

• The risk difference is on the absolute scale: -7.54 excess cases per 100 exposed. Its reciprocal, 13.3, is the number needed to expose for one additional case.

Every measure agrees: smokers died less often. A report written at this point would be internally consistent, statistically significant, and false.

Age

Smoking was less common among the oldest women in the survey, and the oldest women were the ones most likely to die within twenty years. Age is associated with the exposure and independently predicts the outcome, which is the definition of a confounder.

smoking <- epi_strata(
  c(15, 270,  12, 327),
  c(80, 167,  53, 147),
  c(44,   6, 165,  28),
  labels   = c("18-44", "45-64", "65+"),
  exposure = c("Smoker", "Non-smoker"),
  outcome  = c("Dead", "Alive")
)
smoking
#> 18-44
#>               Dead  Alive  Total
#>   Smoker        15    270    285
#>   Non-smoker    12    327    339
#>   Total         27    597    624
#> 
#> 45-64
#>               Dead  Alive  Total
#>   Smoker        80    167    247
#>   Non-smoker    53    147    200
#>   Total        133    314    447
#> 
#> 65+
#>               Dead  Alive  Total
#>   Smoker        44      6     50
#>   Non-smoker   165     28    193
#>   Total        209     34    243

Look at the odds ratio inside each age group before pooling anything:

round(mh_odds_ratio(smoking)$stratum_estimates, 3)
#> 18-44 45-64   65+ 
#> 1.514 1.329 1.244

All three are above 1. The crude estimate was 0.68. Adjustment here does not shift the estimate — it reverses it.

Pooling

options(epibyhand.verbose = 2)
mh_odds_ratio(smoking)

Mantel-Haenszel odds ratio

18-44

45-64

65+

Step 1. Odds ratio within each stratum

\[ \begin{aligned} OR_i &= \frac{a_i d_i}{b_i c_i} \end{aligned} \]

Look at these before pooling. If they disagree substantially the stratifying variable is an effect modifier, and a single pooled number hides the finding rather than reporting it. Test with homogeneity().

Step 2. Stratum contributions

\[ \begin{aligned} & R_i = \frac{a_i d_i}{n_i} \qquad S_i = \frac{b_i c_i}{n_i} \end{aligned} \]

S_i is the weight this stratum carries in the pooled estimate. It is largest where the stratum has the most information, so small or unbalanced strata contribute little.

Step 3. Pooled odds ratio

\[ \begin{aligned} OR_{MH} &= \frac{\sum_i R_i}{\sum_i S_i} \\ &= 39.2393 / 29.0673 \\ &= 1.3499 \end{aligned} \]

Equivalently sum(S_i * OR_i) / sum(S_i): a weighted average of the stratum odds ratios with weights S_i. The pooled value must fall between the smallest and largest stratum estimate; if yours does not, the arithmetic is wrong.

Step 4. Standard error of log(OR_MH), Robins-Breslow-Greenland

\[ \begin{aligned} SE &= \sqrt{\frac{\sum_i P_i R_i}{2R^2} + \frac{\sum_i (P_i S_i + Q_i R_i)}{2RS} + \frac{\sum_i Q_i S_i}{2S^2}} \\ &= \operatorname{sqrt}(0.0062 + 0.015 + 0.0089) \\ &= 0.1734 \end{aligned} \]

P_i = (a_i+d_i)/n_i and Q_i = (b_i+c_i)/n_i are the concordant and discordant proportions in each stratum.

Step 5. 95% confidence interval, lower limit

\[ \begin{aligned} \text{lower} &= \exp\left(\log OR_{MH} - z \cdot SE\right) \\ &= \exp(0.3001 - 1.96 \times 0.1734) \\ &= 0.961 \end{aligned} \]

Step 6. 95% confidence interval, upper limit

\[ \begin{aligned} \text{upper} &= \exp\left(\log OR_{MH} + z \cdot SE\right) \\ &= \exp(0.3001 + 1.96 \times 0.1734) \\ &= 1.8963 \end{aligned} \]

Step 7. Crude odds ratio, ignoring the strata

\[ \begin{aligned} OR_{\text{crude}} &= \frac{A D}{B C} \quad \text{on the collapsed table} \\ &= (139 \times 502) / (443 \times 230) \\ &= 0.6848 \end{aligned} \]

The crude estimate differs from the adjusted one by -49.3%. A change beyond about 10% is the usual working signal that the stratifying variable confounds the association. This is a judgement about the data, not a hypothesis test – do not decide it with a p-value.

Result: \(OR_{MH} = 1.3499\) (95% CI 0.961 to 1.8963)

• Pooling assumes one common odds ratio underlies every stratum. Check that with homogeneity() before reporting this number.

The weights are the part worth dwelling on. S_i is what each stratum contributes, and the pooled estimate is their weighted average — which is why OR_MH must land between the smallest and largest stratum estimate. It does, at 1.35 between 1.24 and 1.51. The crude estimate of 0.68 does not, and could not, because it is not an average of these numbers at all. It is a different quantity that happens to be computed from the same table.

The last step prints the crude estimate beside the adjusted one so the comparison is arithmetic rather than assertion.

Was pooling legitimate?

A single pooled odds ratio only means something if one odds ratio underlies every stratum. If the strata genuinely differ, the stratifying variable is an effect modifier and pooling destroys the finding rather than reporting it.

options(epibyhand.verbose = 1)
homogeneity(smoking)

Breslow-Day test of homogeneity

Step 1. Expected exposed cases in each stratum under a common OR

\[ \begin{aligned} A_i &= \text{root of } (1-\psi)A^2 + \left(N - n_1 - m_1 + \psi(n_1 + m_1)\right)A - \psi n_1 m_1 = 0 \end{aligned} \]

A_i is what cell a would be if this stratum had exactly the pooled odds ratio, holding its margins fixed. Solving a quadratic is the one step here you would not do by hand.

Step 2. Contribution of each stratum to the statistic

\[ \begin{aligned} X^2_i &= \frac{(a_i - A_i)^2}{\operatorname{Var}(A_i)} \end{aligned} \]

A single large contribution means one stratum is driving the result.

Step 3. Statistic, with Tarone’s correction

\[ \begin{aligned} X^2 &= \sum_i X^2_i - \frac{\left(\sum_i a_i - \sum_i A_i\right)^2}{\sum_i \operatorname{Var}(A_i)} \\ &= 0.1182 \end{aligned} \]

Step 4. Reference distribution

\[ \begin{aligned} p &= P\left(\chi^2_{K-1} > X^2\right) \\ &= 0.9426 \end{aligned} \]

Result: \(X^2 = 0.1182\)

• p = 0.9426 on 2 degrees of freedom.
• A large p-value is not evidence that the odds ratios are equal. This test has poor power with small strata, so it will usually fail to reject whether or not effect modification is present. Inspect the stratum-specific estimates as well; they are the more informative thing.

The statistic is small and the p-value large, and the per-stratum contributions show no single stratum straining against the others. Here the assumption is comfortable.

Note the caveat the function prints anyway. A large p-value is not evidence of homogeneity: this test has poor power, and with small strata it will fail to reject almost regardless of the truth. The stratum-specific estimates you looked at earlier remain the more informative thing.

Attributable fractions

Attributable fractions inherit whatever confounding is present, so computing one on the crude table would propagate the error rather than fix it. Work inside a stratum instead, where age is held fixed by construction:

middle <- epi2x2(80, 167, 53, 147,
                 exposure = c("Smoker", "Non-smoker"),
                 outcome  = c("Dead", "Alive"))

options(epibyhand.verbose = 2)
attributable_fraction(middle, among = "population")

Population attributable fraction

Step 1. Risk of Dead among the exposed (Smoker)

\[ \begin{aligned} R_1 &= \frac{a}{a + b} \\ &= 80 / (80 + 167) \\ &= 0.3239 \end{aligned} \]

Step 2. Risk of Dead among the unexposed (Non-smoker)

\[ \begin{aligned} R_0 &= \frac{c}{c + d} \\ &= 53 / (53 + 147) \\ &= 0.265 \end{aligned} \]

Step 3. Risk ratio

\[ \begin{aligned} RR &= \frac{R_1}{R_0} \\ &= 0.3239 / 0.265 \\ &= 1.2222 \end{aligned} \]

Step 4. Risk in the whole population

\[ \begin{aligned} R_t &= \frac{a + c}{n} \\ &= (80 + 53) / 447 \\ &= 0.2975 \end{aligned} \]

Step 5. Population attributable fraction, directly from risks

\[ \begin{aligned} PAF &= \frac{R_t - R_0}{R_t} \\ &= (0.2975 - 0.265) / 0.2975 \\ &= 0.1094 \end{aligned} \]

The share of the population’s risk that would disappear if everyone had the risk of the unexposed.

Step 6. Levin’s formula, from exposure prevalence

\[ \begin{aligned} PAF &= \frac{p(RR - 1)}{1 + p(RR - 1)} \\ &= 0.5526(1.2222 - 1) / (1 + 0.5526(1.2222 - 1)) \\ &= 0.1094 \end{aligned} \]

p = (a + b)/n = 0.5526 is the proportion exposed. Use this form when you have the risk ratio from one study and exposure prevalence from another.

Step 7. Miettinen’s formula, from the proportion of cases exposed

\[ \begin{aligned} PAF &= p_c \cdot \frac{RR - 1}{RR} \\ &= 0.6015 \times (1.2222 - 1) / 1.2222 \\ &= 0.1094 \end{aligned} \]

pc = a/(a + c) = 0.6015 is the proportion of cases who were exposed. All three lines above are the same number. They are one quantity written three ways, not three estimators to choose between.

Step 8. 95% confidence interval

\[ \begin{aligned} \text{CI} &= \text{Levin's formula applied to each limit of the } RR \text{ interval} \end{aligned} \]

This holds exposure prevalence fixed at its observed value, so the interval is slightly too narrow. It is the version you can compute by hand; a delta-method interval that propagates uncertainty in p as well is wider.

Result: \(PAF = 0.1094\) (95% CI -0.0511 to 0.2607)

• Read as: 10.9% of all cases in this population are attributable to the exposure, if the association is causal.
• PAF depends on how common the exposure is, so it does not transfer between populations the way a risk ratio does. A strong risk factor that is rare has a small PAF; a weak one that is universal can have a large one. This is also why PAF from a case-control study needs the exposure prevalence of the source population, not of the controls.

Three formulas for the population attributable fraction circulate in textbooks and students are rarely shown that they are the same quantity. The derivation computes all three — directly from risks, by Levin’s formula from exposure prevalence, and by Miettinen’s from the proportion of cases exposed — and they agree to the last digit, because they are one number written three ways.

Checking a hand calculation

Give check_work() a value and it compares it against the final estimate. When that does not match, it searches every intermediate step for one that does:

options(epibyhand.verbose = 0)
d <- odds_ratio(crude)

check_work(d, 0.3137)
#> Not a match. You gave 0.3137; the final estimate (OR) is 0.6848.
#> 
#> Your value does match step 1: Odds of Dead among the exposed (Smoker).
#>   odds1 = a / b
#> You may have stopped early. The next step is: Odds of Dead among the unexposed (Non-smoker).

The student did not fail; they computed the odds among the exposed and stopped. That is a different problem from an arithmetic slip, and it needs a different sentence from the person teaching them.

Building problem sets

steps_table() returns the whole derivation as a data frame, which is what you want when generating answer keys or rendering the working somewhere this package does not reach:

steps_table(mh_odds_ratio(smoking))[, c("symbol", "label", "result")]
#>     symbol                                                  label    result
#> 1     OR_i                         Odds ratio within each stratum        NA
#> 2     <NA>                                  Stratum contributions        NA
#> 3    OR_MH                                      Pooled odds ratio 1.3499462
#> 4       SE Standard error of log(OR_MH), Robins-Breslow-Greenland 0.1733784
#> 5    lower                   95% confidence interval, lower limit 0.9610291
#> 6    upper                   95% confidence interval, upper limit 1.8962535
#> 7 OR_crude                  Crude odds ratio, ignoring the strata 0.6848366

Output detail is controlled globally with options(epibyhand.verbose = ), where 0 prints the result alone, 1 adds the symbolic formulas, and 2 shows the full worked solution. Set epibyhand.digits to change rounding.

Predictive values

Adjustment is not the only place where a number depends on the population rather than the test. A screening example makes the point sharply:

options(epibyhand.verbose = 1)

screening <- epi2x2(90, 180, 10, 1720,
                    exposure = c("Test positive", "Test negative"),
                    outcome  = c("Diseased", "Healthy"))

predictive_value(screening)

Positive predictive value

Step 1. Sensitivity: of those who are Diseased, the proportion the test finds

\[ \begin{aligned} \text{Sens} &= \frac{TP}{TP + FN} \\ &= 0.9 \end{aligned} \]

Step 2. Specificity: of those who are Healthy, the proportion the test clears

\[ \begin{aligned} \text{Spec} &= \frac{TN}{FP + TN} \\ &= 0.9053 \end{aligned} \]

Step 3. Prevalence in this table

\[ \begin{aligned} \text{Prev} &= \frac{TP + FN}{n} \\ &= 0.05 \end{aligned} \]

Sensitivity and specificity are properties of the test. Everything below also depends on this number.

Step 4. Positive predictive value, by Bayes’ theorem

\[ \begin{aligned} \text{PPV} &= \frac{\text{Sens} \cdot \text{Prev}}{\text{Sens} \cdot \text{Prev} + (1 - \text{Spec})(1 - \text{Prev})} \\ &= 0.3333 \end{aligned} \]

Equivalently TP / (TP + FP) = 90 / (90 + 180) = 0.3333. The two agree because the prevalence used is the one in this table.

Step 5. Standard error of the proportion

\[ \begin{aligned} SE &= \sqrt{\frac{p(1 - p)}{n}} \\ &= 0.0287 \end{aligned} \]

A predictive value is a proportion of those the test called, so n here is the number the test called, not the whole table.

Step 6. 95% confidence interval

\[ \begin{aligned} \text{CI} &= p \pm z \cdot SE \end{aligned} \]

Normal approximation, clamped to [0, 1]. It is the version you can compute by hand; near 0 or 1 an exact or Wilson interval is better behaved.

Result: \(\text{PPV} = 0.3333\) (95% CI 0.2771 to 0.3896)

• Read as: of everyone this test calls positive, 33.3% really are diseased.
• Predictive values are not properties of the test. Move the same test to a population with different prevalence and this number changes, while sensitivity and specificity do not. That is why a test that performs well in a clinic can be useless as a population screen.

A test that is 90% sensitive and 90% specific, applied where 5% of people are actually ill, is wrong two times out of three when it says you are sick.

Sensitivity and specificity are properties of the test. The predictive value is a property of the test and the population, which is why prevalence can be supplied to ask what the same test would do somewhere else:

options(epibyhand.verbose = 0)
vapply(c(0.01, 0.05, 0.20, 0.50),
       function(p) estimate(predictive_value(screening, prevalence = p)),
       numeric(1))
#> [1] 0.0875576 0.3333333 0.7037037 0.9047619

Nothing about the test changed across those four numbers.

In your own documents

Every derivation in this vignette has been typeset as mathematics, with fractions, square roots and subscripts, and none of the code above asks for that. When a derivation is the value of a chunk in a Quarto or R Markdown document, it is rendered this way automatically, to HTML, PDF or Word.

That makes an answer key a matter of writing the questions: the worked solutions come from the same calls a student would make.

For problem sets written in LaTeX directly, to_latex() returns the same derivation as a fragment that needs only \usepackage{amsmath}:

cat(to_latex(odds_ratio(crude), verbose = 1))
#> \medskip\noindent\textbf{Odds ratio}\par\smallskip
#> \noindent\textbf{Step 1.} Odds of Dead among the exposed (Smoker)
#> \[
#> \begin{aligned}
#> \text{odds}_1 &= \frac{a}{b} \\
#> &= 0.3138
#> \end{aligned}
#> \]
#> \medskip
#> 
#> \noindent\textbf{Step 2.} Odds of Dead among the unexposed (Non-smoker)
#> \[
#> \begin{aligned}
#> \text{odds}_0 &= \frac{c}{d} \\
#> &= 0.4582
#> \end{aligned}
#> \]
#> \medskip
#> 
#> \noindent\textbf{Step 3.} Odds ratio
#> \[
#> \begin{aligned}
#> OR &= \frac{\text{odds}_1}{\text{odds}_0} \\
#> &= 0.6848
#> \end{aligned}
#> \]
#> {\small\emph{Equivalently OR = ad / bc, which is why the odds ratio is the same whether you condition on exposure or on outcome. That symmetry is what makes it usable in case-control studies.}}\par
#> \medskip
#> 
#> \noindent\textbf{Step 4.} Standard error of log(OR), Woolf's method
#> \[
#> \begin{aligned}
#> SE &= \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \\
#> &= 0.1257
#> \end{aligned}
#> \]
#> \medskip
#> 
#> \noindent\textbf{Step 5.} 95\% confidence interval, lower limit
#> \[
#> \begin{aligned}
#> \text{lower} &= \exp\left(\log OR - z \cdot SE\right) \\
#> &= 0.5353
#> \end{aligned}
#> \]
#> \medskip
#> 
#> \noindent\textbf{Step 6.} 95\% confidence interval, upper limit
#> \[
#> \begin{aligned}
#> \text{upper} &= \exp\left(\log OR + z \cdot SE\right) \\
#> &= 0.8761
#> \end{aligned}
#> \]
#> \medskip
#> 
#> \noindent\textbf{Result:} $OR = 0.6848$ (95\% CI 0.5353 to 0.8761)\par
#> \begin{itemize}
#> \item These data are tabulated as a cohort, where risk is estimable. The risk ratio here is 0.7601 against an odds ratio of 0.6848 -- the odds ratio is the more extreme of the two, and always will be. The two converge only when the outcome is rare; baseline risk here is 0.3142.
#> \end{itemize}

standalone = TRUE wraps it in a complete document. If you would rather keep the console style in a document, set options(epibyhand.knit = "console").

Scope

This package covers methods a student can compute by hand on paper. The boundary is deliberate. It is what keeps the package small enough to stay correct without constant maintenance, and it is why there is no regression modelling here — once the estimate comes from an iterative fit there is no hand calculation left to check, and a printed “derivation” would be decoration rather than instruction.