---
title: "Working through a confounded analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Working through a confounded analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(epibyhand)
```

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.

## 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.

```{r}
crude <- epi2x2(139, 443, 230, 502,
                exposure = c("Smoker", "Non-smoker"),
                outcome  = c("Dead", "Alive"))
crude
```

## The crude analysis

```{r}
odds_ratio(crude)
```

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:

```{r}
options(epibyhand.verbose = 1)
risk_ratio(crude)
risk_difference(crude)
```

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.

```{r}
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
```

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

```{r}
round(mh_odds_ratio(smoking)$stratum_estimates, 3)
```

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

## Pooling

```{r}
options(epibyhand.verbose = 2)
mh_odds_ratio(smoking)
```

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.

```{r}
options(epibyhand.verbose = 1)
homogeneity(smoking)
```

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:

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

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

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:

```{r}
options(epibyhand.verbose = 0)
d <- odds_ratio(crude)

check_work(d, 0.3137)
```

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:

```{r}
steps_table(mh_odds_ratio(smoking))[, c("symbol", "label", "result")]
```

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.

```{r, include = FALSE}
options(epibyhand.verbose = 2, epibyhand.digits = 4)
```

## 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.
