The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Statistical software gives you the answer. In a methods course the
answer is the least interesting part of the calculation: a student who
can produce 2.14 without being able to say where it came
from has learned nothing that will survive the exam.
Every function in this package returns the answer together with the reasoning that produced it. This vignette works one dataset the whole way through, in roughly the order a first course would, and the derivations are the point rather than an aside.
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.
odds_ratio(crude)
#> ── Odds ratio ──────────────────────────────────────────────────────────────
#>
#> Data
#> Dead Alive Total
#> Smoker 139 443 582
#> Non-smoker 230 502 732
#> Total 369 945 1314
#>
#> Step 1 Odds of Dead among the smoker
#> odds1 = a / b
#> = 139 / 443
#> = 0.3138
#>
#> Step 2 Odds of Dead among the non-smoker
#> odds0 = c / d
#> = 230 / 502
#> = 0.4582
#>
#> Step 3 Odds ratio
#> OR = odds1 / odds0
#> = 0.3138 / 0.4582
#> = 0.6848
#> # 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
#> SE = sqrt(1/a + 1/b + 1/c + 1/d)
#> = sqrt(1/139 + 1/443 + 1/230 + 1/502)
#> = 0.1257
#>
#> Step 5 95% confidence interval, lower limit
#> lower = exp(log(OR) - z * SE)
#> = exp(-0.3786 - 1.96 * 0.1257)
#> = 0.5353
#>
#> Step 6 95% confidence interval, upper limit
#> upper = exp(log(OR) + z * SE)
#> = exp(-0.3786 + 1.96 * 0.1257)
#> = 0.8761
#>
#>
#> Result
#> OR = 0.6848 (95% CI 0.5353 to 0.8761)
#>
#> Notes
#> 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 smoker
#> R1 = a / (a + b)
#> = 0.2388
#>
#> Step 2 Risk of Dead among the non-smoker
#> R0 = c / (c + d)
#> = 0.3142
#>
#> Step 3 Risk ratio
#> RR = R1 / R0
#> = 0.7601
#>
#> Step 4 Standard error of log(RR)
#> SE = sqrt(1/a - 1/(a+b) + 1/c - 1/(c+d))
#> = 0.092
#> # 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
#> lower = exp(log(RR) - z * SE)
#> = 0.6347
#>
#> Step 6 95% confidence interval, upper limit
#> upper = exp(log(RR) + z * SE)
#> = 0.9102
#>
#>
#> Result
#> RR = 0.7601 (95% CI 0.6347 to 0.9102)
#>
#> Notes
#> 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 smoker
#> R1 = a / (a + b)
#> = 0.2388
#>
#> Step 2 Risk of Dead among the non-smoker
#> R0 = c / (c + d)
#> = 0.3142
#>
#> Step 3 Risk difference
#> RD = R1 - R0
#> = -0.0754
#>
#> Step 4 Standard error of RD
#> SE = sqrt(R1(1-R1)/(a+b) + R0(1-R0)/(c+d))
#> = 0.0246
#> # 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
#> lower = RD - z * SE
#> = -0.1237
#>
#> Step 6 95% confidence interval, upper limit
#> upper = RD + z * SE
#> = -0.0271
#>
#>
#> Result
#> RD = -0.0754 (95% CI -0.1237 to -0.0271)
#>
#> Notes
#> 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.
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 243Look at the odds ratio inside each age group before pooling anything:
All three are above 1. The crude estimate was 0.68. Adjustment here does not shift the estimate — it reverses it.
options(epibyhand.verbose = 2)
mh_odds_ratio(smoking)
#> ── Mantel-Haenszel odds ratio ──────────────────────────────────────────────
#>
#> Data
#> 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
#>
#> Step 1 Odds ratio within each stratum
#> OR_i = (a_i * d_i) / (b_i * c_i)
#>
#> Stratum a b c d n_i OR_i
#> ------- -- --- --- --- --- ------
#> 18-44 15 270 12 327 624 1.5139
#> 45-64 80 167 53 147 447 1.3287
#> 65+ 44 6 165 28 243 1.2444
#>
#> # 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
#> R_i = a_i*d_i/n_i S_i = b_i*c_i/n_i
#>
#> Stratum R_i S_i
#> ------- ------- -------
#> 18-44 7.8606 5.1923
#> 45-64 26.3087 19.8009
#> 65+ 5.07 4.0741
#>
#> # 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
#> OR_MH = sum(R_i) / sum(S_i)
#> = 39.2393 / 29.0673
#> = 1.3499
#> # 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
#> SE = sqrt( sum(P_i R_i)/(2R^2) + sum(P_i S_i + Q_i R_i)/(2RS) + sum(Q_i S_i)/(2S^2) )
#> = sqrt(0.0062 + 0.015 + 0.0089)
#> = 0.1734
#>
#> Stratum P_i Q_i
#> ------- ------ ------
#> 18-44 0.5481 0.4519
#> 45-64 0.5078 0.4922
#> 65+ 0.2963 0.7037
#>
#> # 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
#> lower = exp(log(OR_MH) - z * SE)
#> = exp(0.3001 - 1.96 * 0.1734)
#> = 0.961
#>
#> Step 6 95% confidence interval, upper limit
#> upper = exp(log(OR_MH) + z * SE)
#> = exp(0.3001 + 1.96 * 0.1734)
#> = 1.8963
#>
#> Step 7 Crude odds ratio, ignoring the strata
#> OR_crude = (A * D) / (B * C) on the collapsed table
#> = (139 * 502) / (443 * 230)
#> = 0.6848
#> # 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)
#>
#> Notes
#> 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.
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
#> A_i = root of (1-psi)A^2 + (N - n1 - m1 + psi(n1+m1))A - psi*n1*m1 = 0
#> # 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
#> X2_i = (a_i - A_i)^2 / Var(A_i)
#> # A single large contribution means one stratum is driving the
#> # result.
#>
#> Step 3 Statistic, with Tarone's correction
#> X2 = sum(X2_i) - (sum(a_i) - sum(A_i))^2 / sum(Var(A_i))
#> = 0.1182
#>
#> Step 4 Reference distribution
#> p = P(chi-squared with K - 1 df > X2)
#> = 0.9426
#>
#>
#> Result
#> X2 = 0.1182
#>
#> Notes
#> 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 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 ────────────────────────────────────────
#>
#> Data
#> Dead Alive Total
#> Smoker 80 167 247
#> Non-smoker 53 147 200
#> Total 133 314 447
#>
#> Step 1 Risk of Dead among the smoker
#> R1 = a / (a + b)
#> = 80 / (80 + 167)
#> = 0.3239
#>
#> Step 2 Risk of Dead among the non-smoker
#> R0 = c / (c + d)
#> = 53 / (53 + 147)
#> = 0.265
#>
#> Step 3 Risk ratio
#> RR = R1 / R0
#> = 0.3239 / 0.265
#> = 1.2222
#>
#> Step 4 Risk in the whole population
#> Rt = (a + c) / n
#> = (80 + 53) / 447
#> = 0.2975
#>
#> Step 5 Population attributable fraction, directly from risks
#> PAF = (Rt - R0) / Rt
#> = (0.2975 - 0.265) / 0.2975
#> = 0.1094
#> # 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
#> PAF = p(RR - 1) / (1 + p(RR - 1))
#> = 0.5526(1.2222 - 1) / (1 + 0.5526(1.2222 - 1))
#> = 0.1094
#> # 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
#> PAF = pc * (RR - 1) / RR
#> = 0.6015 * (1.2222 - 1) / 1.2222
#> = 0.1094
#> # 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
#> CI = Levin's formula applied to each limit of the RR interval
#>
#> limit RR PAF
#> ----- ----- -------
#> lower 0.912 -0.0511
#> upper 1.638 0.2607
#>
#> # 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)
#>
#> Notes
#> 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.
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 smoker.
#> odds1 = a / b
#> You may have stopped early. The next step is: Odds of Dead among the 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.
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.6848366Output 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.
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.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.