{
  "format": "srwb",
  "version": "1.0",
  "notebook": {
    "name": "Simpson's Paradox",
    "cells": [
      {
        "type": "markdown",
        "language": "markdown",
        "name": "intro",
        "code": "# Simpson's Paradox: When Every Group Disagrees with the Total\n\n> **All data in this workbook is synthetic/fictional**, constructed for classroom illustration. No real patients, products, or treatments are involved.\n\n## Learning objectives\nBy the end of this notebook you should be able to:\n\n1. Construct a small stratified 2×2 dataset that exhibits Simpson's paradox.\n2. Compute success proportions and **relative risks (RR = p_A / p_B)** within each stratum and after pooling.\n3. Explain **why** the aggregate comparison can reverse the within-stratum comparison (confounding by stratum size / a lurking variable).\n4. Verify the same numbers independently in a second language (R and JavaScript).\n\n## Prediction question (answer before running the computation cell)\nSuppose Treatment A beats Treatment B **in every stratum** (higher success rate, RR > 1 everywhere). Is it possible that when you ignore the strata and pool all the data, Treatment A comes out **worse** than Treatment B?\n\n- Write down your guess: **Yes / No**, and one sentence of reasoning.\n- If you answered *Yes*: what must be true about the **stratum sizes** of A and B for the reversal to appear?\n- Try computing the pooled rates by hand once you see the counts in the next cell.\n\n## Runtime and run order\n\nRun the named cells in order: `data` (R), `compute` (R), `viz` (R), then `js_check` (JavaScript). The R cells share runtime variables; showing a saved cell is not the same as running it. This catalog copy contains no saved outputs.\n\nThe R example uses base R only: no package installation is needed, but SciREPL may need to download or initialize its webR runtime the first time. JavaScript runs in the browser. The JavaScript cell independently recalculates the original fixed counts; it does not read the R session. If you change the dataset, update both languages' input counts and their checks before comparing them."
      },
      {
        "type": "code",
        "language": "r",
        "name": "data",
        "code": "# ---- SYNTHETIC / FICTIONAL DATA (not real) ----\n# Two strata; two treatments A and B; binary outcome (success/failure).\n# Constructed so that A > B within each stratum, but B > A when pooled.\n# Mechanism: A treats mostly Severe patients (high-risk), B treats mostly Mild patients.\n\ndat <- data.frame(\n  stratum   = rep(c(\"Mild\", \"Severe\"), each = 4),\n  treatment = rep(rep(c(\"A\", \"B\"), each = 2), 2),\n  outcome   = rep(c(\"success\", \"failure\"), 4),\n  count = c(\n    # Mild stratum\n    85, 15,       # A: 85 success, 15 failure   -> p = 0.85  (n = 100)\n    350, 150,     # B: 350 success, 150 failure -> p = 0.70  (n = 500)\n    # Severe stratum\n    30, 470,      # A: 30 success, 470 failure  -> p = 0.06  (n = 500)\n    2,  98        # B:  2 success, 98 failure   -> p = 0.02  (n = 100)\n  )\n)\ndat"
      },
      {
        "type": "markdown",
        "language": "markdown",
        "name": "pause",
        "code": "## ⏸️ Pause — predict and calculate\n\nThe dataset above is now defined. **Do not run the next cell yet.**\n\n1. From the counts, compute by hand the success proportion of A and B in the *Mild* stratum, and in the *Severe* stratum. Which treatment wins in each stratum?\n2. Now pool: add up **all** successes and **all** failures across both strata for each treatment, and compute the pooled proportions. Which treatment wins overall?\n3. Compute the three relative risks RR = p_A / p_B (Mild, Severe, pooled). Do they all point the same way?\n4. Before revealing the answers: *why* does A look better everywhere yet worse overall? (Hint: which treatment treats mostly **Severe** patients, and what does that do to its overall average?)\n\nWrite your predictions down, then run the next cell to check."
      },
      {
        "type": "code",
        "language": "r",
        "name": "compute",
        "code": "# ---- Computation from the dataset (nothing hard-coded below) ----\n\n# Derive strata from the data cell: run data before this cell.\nstrata <- unique(dat$stratum)\n\nget_p <- function(s, t) {\n  n_ok <- dat$count[dat$stratum == s & dat$treatment == t & dat$outcome == \"success\"]\n  n_no <- dat$count[dat$stratum == s & dat$treatment == t & dat$outcome == \"failure\"]\n  c(success = n_ok, failure = n_no, n = n_ok + n_no, p = n_ok / (n_ok + n_no))\n}\n\nrr <- function(pA, pB) pA / pB\n\n# --- Within-stratum results ---\nfor (s in strata) {\n  a <- get_p(s, \"A\"); b <- get_p(s, \"B\")\n  cat(sprintf(\"--- Stratum: %s ---\\n\", s))\n  cat(sprintf(\"  A: success = %d, failure = %d  (p_A = %.4f, n = %d)\\n\", a[\"success\"], a[\"failure\"], a[\"p\"], a[\"n\"]))\n  cat(sprintf(\"  B: success = %d, failure = %d  (p_B = %.4f, n = %d)\\n\", b[\"success\"], b[\"failure\"], b[\"p\"], b[\"n\"]))\n  cat(sprintf(\"  RR(A/B) in %s = %.4f\\n\\n\", s, rr(a[\"p\"], b[\"p\"])))\n}\n\n# --- Pooled results (collapse over strata) ---\nagg <- aggregate(count ~ treatment + outcome, data = dat, FUN = sum)\npooled <- function(t) {\n  s <- agg$count[agg$treatment == t & agg$outcome == \"success\"]\n  f <- agg$count[agg$treatment == t & agg$outcome == \"failure\"]\n  c(success = s, failure = f, n = s + f, p = s / (s + f))\n}\npA <- pooled(\"A\"); pB <- pooled(\"B\")\ncat(\"--- Pooled (strata collapsed) ---\\n\")\ncat(sprintf(\"  A: success = %d, failure = %d  (p_A = %.4f, n = %d)\\n\", pA[\"success\"], pA[\"failure\"], pA[\"p\"], pA[\"n\"]))\ncat(sprintf(\"  B: success = %d, failure = %d  (p_B = %.4f, n = %d)\\n\", pB[\"success\"], pB[\"failure\"], pB[\"p\"], pB[\"n\"]))\ncat(sprintf(\"  Pooled RR(A/B) = %.4f\\n\\n\", rr(pA[\"p\"], pB[\"p\"])))\n\n# --- Proofs of the paradox ---\nrr_mild    <- rr(get_p(\"Mild\",   \"A\")[\"p\"], get_p(\"Mild\",   \"B\")[\"p\"])\nrr_severe  <- rr(get_p(\"Severe\", \"A\")[\"p\"], get_p(\"Severe\", \"B\")[\"p\"])\nrr_pooled  <- rr(pA[\"p\"], pB[\"p\"])\n\nstopifnot(\n  rr_mild   > 1,   # A beats B within Mild\n  rr_severe > 1,   # A beats B within Severe\n  rr_pooled < 1    # ... yet A loses to B overall -> Simpson's paradox\n)\ncat(\"✔ All stopifnot checks passed: RR_mild > 1, RR_severe > 1, but RR_pooled < 1.\")"
      },
      {
        "type": "code",
        "language": "r",
        "name": "viz",
        "code": "# ---- Visualization: success rates by stratum and pooled (base R only) ----\n\nlabels <- c(\"Mild\", \"Severe\", \"Pooled\")\npA_v <- c(get_p(\"Mild\", \"A\")[\"p\"], get_p(\"Severe\", \"A\")[\"p\"], pA[\"p\"])\npB_v <- c(get_p(\"Mild\", \"B\")[\"p\"], get_p(\"Severe\", \"B\")[\"p\"], pB[\"p\"])\n\nm <- rbind(A = pA_v, B = pB_v)\ncolnames(m) <- labels\nop <- par(mar = c(4.1, 4.1, 4.1, 1.1))\nbarplot(m, beside = TRUE, ylim = c(0, 0.9),\n        col = c(\"#4472C4\", \"#ED7D31\"),\n        main = paste(strwrap(\"Simpson's Paradox: success rate by group (synthetic data)\", width = 42), collapse = \"\\n\"), cex.main = 1.0,\n        ylab = \"Success proportion\", legend.text = c(\"Treatment A\", \"Treatment B\"),\n        args.legend = list(x = \"topright\", bty = \"n\"))\nabline(h = 0, col = \"grey40\")\npar(op)"
      },
      {
        "type": "code",
        "language": "javascript",
        "name": "js_check",
        "code": "// Independent recalculation from the underlying counts (SYNTHETIC data)\nconst counts = {\n  Mild:   { A: { s: 85,  f: 15 },  B: { s: 350, f: 150 } },\n  Severe: { A: { s: 30,  f: 470 }, B: { s: 2,   f: 98 } }\n};\n\nconst res = {};\nfor (const st of Object.keys(counts)) {\n  const g = counts[st];\n  for (const t of [\"A\", \"B\"]) g[t].p = g[t].s / (g[t].s + g[t].f);\n  res[st] = { A: g.A, B: g.B, RR: g.A.p / g.B.p };\n}\n\n// Pool by summing successes and failures across strata\nconst pooled = { A: { s: 0, f: 0 }, B: { s: 0, f: 0 } };\nfor (const st of Object.keys(counts)) {\n  pooled.A.s += counts[st].A.s; pooled.A.f += counts[st].A.f;\n  pooled.B.s += counts[st].B.s; pooled.B.f += counts[st].B.f;\n}\nfor (const t of [\"A\", \"B\"]) pooled[t].p = pooled[t].s / (pooled[t].s + pooled[t].f);\nres.Pooled = { A: pooled.A, B: pooled.B, RR: pooled.A.p / pooled.B.p };\n\n// Readable summary\nconsole.log(\"Independent JavaScript recalculation (synthetic data):\");\nfor (const k of Object.keys(res)) {\n  const r = res[k];\n  console.log(\n    k.padEnd(7) +\n    `  p_A = ${r.A.p.toFixed(4)} (${r.A.s}/${r.A.s + r.A.f})` +\n    `  p_B = ${r.B.p.toFixed(4)} (${r.B.s}/${r.B.s + r.B.f})` +\n    `  RR = ${r.RR.toFixed(4)}`\n  );\n}\n\n// Check the original fixed counts independently; this does not read R state.\nconst expect = {\n  Mild:   { pA: 85/100,  pB: 350/500 },\n  Severe: { pA: 30/500,  pB: 2/100 },\n  Pooled: { pA: 115/600, pB: 352/600 }\n};\nconst eps = 1e-12;\nfunction near(a, b) { return Math.abs(a - b) < eps; }\n\nif (!near(res.Mild.A.p,   expect.Mild.pA)   || !near(res.Mild.B.p,   expect.Mild.pB))   throw new Error(\"Mild proportions disagree with the expected input counts\");\nif (!near(res.Severe.A.p, expect.Severe.pA) || !near(res.Severe.B.p, expect.Severe.pB)) throw new Error(\"Severe proportions disagree with the expected input counts\");\nif (!near(res.Pooled.A.p, expect.Pooled.pA) || !near(res.Pooled.B.p, expect.Pooled.pB)) throw new Error(\"Pooled proportions disagree with the expected input counts\");\nif (!(res.Mild.RR > 1 && res.Severe.RR > 1 && res.Pooled.RR < 1))\n  throw new Error(\"Paradox structure not reproduced: expected RR_mild > 1, RR_severe > 1, RR_pooled < 1\");\n\nconsole.log(\"\\n✔ Original-count JavaScript recalculation and paradox checks passed.\");"
      },
      {
        "type": "markdown",
        "language": "markdown",
        "name": "explanation",
        "code": "## Why the reversal happens\n\n**The paradox.** In both strata, A's success rate exceeds B's (RR > 1 in each). Yet when the strata are collapsed, B's overall rate exceeds A's. Both statements are true of the same data — aggregation changed the conclusion.\n\n**The mechanism: confounding by severity.** Severity is a *lurking variable* that shapes both treatment assignment and outcome:\n\n- Severe patients have very low success rates under **either** treatment (6% for A, 2% for B).\n- In this synthetic cohort, **A treats mostly Severe patients** (500 of 600), while **B treats mostly Mild patients** (500 of 600).\n- So A's overall average is dragged down by its heavy Severe case mix, while B's is propped up by its Mild-heavy mix — even though A is better *within every stratum*.\n\nFormally, each treatment's pooled rate is a weighted average of its stratum-specific rates:\n\n$$p_T^{\\,pooled} = w_{T,\\text{Mild}}\\,p_{T,\\text{Mild}} + w_{T,\\text{Severe}}\\,p_{T,\\text{Severe}}, \\qquad w_{T,s} = \\frac{n_{T,s}}{n_T}$$\n\nThe weights $w_{T,s}$ **differ between treatments** (A: Mild weight 1/6 vs B: 5/6), so the pooled difference mixes a within-stratum association with a case-mix difference. When the imbalance in weights is strong enough relative to the gap in stratum rates, the sign flips. Stratifying (or standardizing both treatments to the same stratum weights) removes *this particular measured case-mix distortion*: it shows that **A is better in every stratum of this synthetic data**. It does **not** by itself isolate \"the treatment effect\" or prove A is causally better — a causal conclusion would require an appropriate design (e.g. randomization) and assumptions about confounders *other* than severity that are not represented in this two-variable illustration.\n\n**Why the pooled comparison is misleading here.** The pooled number answers \"which treatment does better *for this particular case mix*?\" — a question dominated by who got assigned where — instead of \"which treatment is better for a patient of a given severity?\" A stratified or common-weight standardized comparison separates this measured case-mix distortion from the within-stratum association. Neither comparison, by itself, establishes a causal treatment effect.\n\n## Three discussion questions\n\n1. **Design:** In this synthetic dataset, how would random assignment reduce the systematic case-mix imbalance? Why might stratified randomization be helpful, and why can chance imbalance still occur?\n2. **Interpretation:** A hospital reports its new therapy has a higher *overall* survival rate than the standard. Before concluding the new therapy is better, which strata (age, baseline severity, comorbidities) would you demand to see?\n3. **Decision:** The pooled RR is < 1 but each stratum RR is > 1. Which estimate is more relevant for a clinician treating one known-severe patient, and why? What additional evidence (study design, other confounders, replication) would be needed before treating that estimate as causal guidance for treatment decisions?"
      }
    ]
  }
}