pacman::p_load(tidyverse)
paths <- c("../lab1/data/derived/analytic_dataset.csv", # Lab 1, completed in place
"data/analytic_dataset.csv") # or the CANVAS copy
path <- paths[file.exists(paths)][1]
if (is.na(path)) {
stop("analytic_dataset.csv not found. Complete Lab 1 first, or download ",
"the dataset from CANVAS into labs/lab2/data/.")
}
lab1 <- read_csv(path, show_col_types = FALSE)
nrow(lab1) # expect 1,292Lab 2 — Three Ways to Sample Controls
EPID 785R · Regression and Study Design
Introduction
In Lab 1 you built an analytic dataset for a target trial emulation: 1,292 ICU patients classified at a 60-minute landmark as early vasopressor initiators or not, followed to hospital discharge, with in-hospital death as the outcome. Because you constructed that cohort yourself, you can enumerate the entire source population. That is, every risk, rate, and odds is directly computable.
This lab uses those data to demonstrate the case-control study. We will draw three different case-control samples from this cohort you created, keeping all the cases each time and changing only how the controls are chosen:
- Cumulative (survivor) sampling are controls drawn from those who never became cases at the end of follow up;
- Case-base sampling are controls drawn from the entire cohort at baseline, blind to what happens later;
- Incidence-density (risk-set) sampling are controls drawn from those still under observation at each case’s event time.
Each sample gets an identical analysis, which is the exposure odds ratio, a cross-product of four cell counts. The point of the lab is that this one estimator, effectively the case-control odds ratio, returns a different cohort parameter under each scheme: the cohort odds ratio, the cohort risk ratio, and the cohort rate ratio, respectively. Because we hold the full cohort, we can check each claim against the truth.
By the end you should be able to:
- compute the 2×2-table risk ratio, odds ratio, and (using person-time) rate ratio directly from a fully enumerated cohort;
- implement all three control-sampling schemes;
- state which cohort parameter each scheme’s exposure odds ratio estimates, and check that claim empirically;
- explain what case-control data alone can’t deliver (absolute risks), and what restores it.
Deliverables (yours to keep — nothing is submitted): this document rendered to HTML with your answers written in the answer blocks, plus the three case-control datasets your code writes to data/ in Step 6 — Week 11’s lab picks those files back up. Once every chunk runs without error, change eval: false to eval: true in the YAML header and re-render, so your HTML shows your results.
1 Before you start (~5 min)
- Open
lab2.Rprojin RStudio. - This lab needs your Lab 1 analytic dataset. If
labs/lab1/sits next tolabs/lab2/and you completed the pipeline, the chunk below finds it automatically. Otherwise, downloadanalytic_dataset.csvfrom CANVAS intolabs/lab2/data/. - Packages: the tidyverse only (
pacman::p_load(tidyverse)installs it if needed).
Two working rules:
- Chunks marked
TODOcontain???placeholders — invalid code, so the chunk errors until you complete it. - Do not change the seeds. Every sampling chunk sets its own seed (
set.seed(...)at the top of the chunk), so your draws match the checkpoints exactly and you can re-run any chunk on its own without disturbing the others. If a checkpoint fails, the bug is in your code, not in the randomness.
2 Step 1 — The cohort answers, computed directly (~12 min)
Before sampling anything, establish the truth we will try to recover. Eighteen stays have a missing outcome (blank hospital discharge status, kept as NA in Lab 1 — all 18 are in the no-early-initiation group). A case-control study needs every sampled person classifiable as case or non-case, so we set them aside here and document the decision.
Complete TODO (1a) — the analytic cohort and its 2×2 counts.
# TODO (1a): keep the stays whose outcome is recorded
coh <- lab1 |> filter(???)
n_coh <- nrow(coh) # expect 1,274
a1 <- sum(coh$early_vasopressor == 1 & coh$hospital_death == 1) # exposed cases
b1 <- sum(coh$early_vasopressor == 1 & coh$hospital_death == 0) # exposed non-cases
a0 <- sum(coh$early_vasopressor == 0 & coh$hospital_death == 1) # unexposed cases
b0 <- sum(coh$early_vasopressor == 0 & coh$hospital_death == 0) # unexposed non-cases
c(n_coh = n_coh, a1 = a1, b1 = b1, a0 = a0, b0 = b0)Checkpoint: 1,274 stays; the 2×2 table is 6 / 21 (exposed) and 130 / 1,117 (unexposed). One Lab 1 lesson worth re-seeing in these counts: the unexposed group contains 107 patients who did start a vasopressor, but after the 60-minute window. “No early initiation” is a strategy, not “never treated.”
Now compute the three cohort parameters that the three sampling schemes will target. Complete TODO (1b) (risks, risk ratio, odds ratio) and TODO (1c) (person-time and the rate ratio; followup_days is each person’s time under observation, so summing it within exposure groups gives person-days at risk).
# TODO (1b): risks and their ratios
R1 <- a1 / ??? # risk in the exposed
R0 <- a0 / ??? # risk in the unexposed
rr_cohort <- ??? / ??? # cohort risk ratio
or_cohort <- (a1 / b1) / (a0 / b0) # cohort odds ratio (provided)
round(c(R1 = R1, R0 = R0, RR = rr_cohort, OR = or_cohort), 3)
# TODO (1c): person-time, rates, rate ratio
pt1 <- sum(coh$followup_days[coh$early_vasopressor == ???])
pt0 <- sum(coh$followup_days[coh$early_vasopressor == ???])
rate1 <- a1 / pt1
rate0 <- a0 / pt0
ratr_cohort <- ??? / ???
round(c(pt1 = pt1, pt0 = pt0,
rate1_per100pd = 100 * rate1, rate0_per100pd = 100 * rate0,
rate_ratio = ratr_cohort), 3)Checkpoints: R1 = 0.222, R0 = 0.104; risk ratio = 2.13; odds ratio = 2.46. Person-days 168.3 (exposed) and 7,310.1 (unexposed); rates 3.56 vs. 1.78 per 100 person-days; rate ratio = 2.00.
These are three different summaries in one cohort, and they differ visibly because the outcome is not rare (10–22% risks). Each sampling scheme below will recover exactly one of them.
Finally, the ingredients every scheme shares: all 136 cases, a target of 4 controls per case, and one helper function.
cases <- coh |> filter(hospital_death == 1) # every scheme keeps all cases
n_cases <- nrow(cases)
n_controls <- 4 * n_cases # 4 controls per case = 544
exp_cases <- sum(cases$early_vasopressor) # 6 exposed cases
# exposure odds ratio of a case series vs. a control series
caco_or <- function(exp_cases, n_cases, exp_controls, n_controls) {
(exp_cases / (n_cases - exp_cases)) /
(exp_controls / (n_controls - exp_controls))
}(Why 4 controls per case, and not 1? With only ~2% exposure prevalence, a 1:1 study would expect just 2–3 exposed controls, which will lead to huge sampling variability. Four per case quadruples that expectation at modest cost.)
Q1. The cohort odds ratio (2.46) sits farther from 1 than the cohort risk ratio (2.13). Why, in this cohort? Under what circumstance would the two be nearly equal?
YOUR ANSWER:
3 Step 2 — Scheme A: cumulative (survivor) sampling (~10 min)
Wait until follow-up is over, then draw controls from the people who never became cases. Complete TODO (2a) — the control pool.
set.seed(153)
# TODO (2a): the pool = cohort members who did NOT experience the outcome
pool_A <- coh |> filter(???)
ctrl_A <- pool_A[sample(nrow(pool_A), n_controls), ]
exp_A <- sum(ctrl_A$early_vasopressor)
or_A <- caco_or(exp_cases, n_cases, exp_A, n_controls)
c(exp_controls_A = exp_A, or_A = round(or_A, 2))Checkpoint: 10 exposed controls; OR_A = 2.46.
Compare that to Step 1: your case-control study of 680 people (136 cases, 544 controls) just reproduced the full cohort’s odds ratio, because a control series drawn from the non-cases stands in for the non-case cells of the cohort 2×2 table, and the sampling fractions cancel out of the cross-product.
Q2. In your Scheme-A sample, 136 of 680 people (20%) are cases, yet the cohort’s risk is 10.7%. Which number is “wrong”? What does the 20% actually reflect, and what family of quantities does this tell you a case-control sample alone can never estimate?
YOUR ANSWER:
4 Step 3 — Scheme B: case-base sampling (~8 min)
Now change one line. Instead of sampling controls from the non-cases at the end of follow-up, draw them from the entire cohort at baseline. Some sampled controls will therefore go on to become cases. It is the design: a baseline sample of the cohort represents the denominators of the risks (everyone who started), not just the survivors.
Complete TODO (3a).
set.seed(353)
# TODO (3a): the pool = the entire cohort at baseline
pool_B <- ???
ctrl_B <- pool_B[sample(nrow(pool_B), n_controls), ]
exp_B <- sum(ctrl_B$early_vasopressor)
fut_B <- sum(ctrl_B$hospital_death == 1) # controls who later become cases
or_B <- caco_or(exp_cases, n_cases, exp_B, n_controls)
c(exp_controls_B = exp_B, future_cases_B = fut_B, or_B = round(or_B, 2))Checkpoint: 11 exposed controls; 61 of the 544 controls later become cases; OR_B = 2.24 — your draw’s estimate of the cohort risk ratio (2.13). This case-control odds ratio estimates a risk ratio not its odds ratio, with no rare-disease approximation anywhere in sight.
Why are they different? Why is the case-control odds ratio 2.24, but the risk ratio 2.13? This is a sampling variability issue. Imbalanced exposure and outcome makes estimating the cohort risk ratio with the case-control odds ratio more difficult. We’ll see more about this below.
Q3. A colleague reviewing your code insists the 61 “contaminated” controls who became cases must be deleted. If you obeyed, what scheme would your design turn into, and toward which cohort parameter would your odds ratio drift? What does keeping them buy?
YOUR ANSWER:
5 Step 4 — Scheme C: incidence-density (risk-set) sampling (~12 min)
The third scheme samples controls longitudinally: at each case’s event time, draw controls from the people still under observation at that moment — here, still in hospital, i.e., followup_days at least as large as the case’s. A person’s chance of ever being sampled is then proportional to the time they spend at risk, so the control series represents person-time, and the exposure odds ratio estimates the rate ratio.
Two features to expect before you run it: the same person may be sampled at several event times, and a sampled control may later be a case — both correct, for the same denominators-are-person-time reason. One quirk of our real data: at the latest death (51.7 days after the landmark) only one other patient is still in hospital, so that risk set yields 1 control, not 4, and the control series totals 541 rather than 544.
Complete TODO (4a) — the risk-set condition.
set.seed(553)
picks <- vector("list", n_cases)
for (j in seq_len(n_cases)) {
t_j <- cases$followup_days[j]
# TODO (4a): at risk at t_j = still under observation at t_j
# (followup_days >= t_j), and not the index case itself
rs <- which(coh$??? >= t_j &
coh$patientunitstayid != cases$patientunitstayid[j])
# sample.int() sidesteps R's sample(x, k) surprise when length(x) == 1
picks[[j]] <- rs[sample.int(length(rs), min(4, length(rs)))]
}
idx_C <- unlist(picks)
ctrl_C <- coh[idx_C, ]
n_C <- nrow(ctrl_C)
exp_C <- sum(ctrl_C$early_vasopressor)
fut_C <- sum(ctrl_C$hospital_death == 1)
rep_C <- sum(duplicated(idx_C))
or_C <- caco_or(exp_cases, n_cases, exp_C, n_C)
c(n_controls_C = n_C, exp_controls_C = exp_C, or_C = round(or_C, 2),
future_cases_C = fut_C, sampled_more_than_once_C = rep_C)Checkpoint: 541 controls; 12 exposed; OR_C = 2.03 — right at the cohort rate ratio (2.00); 57 sampled controls later die; 138 selections are repeats of someone already sampled.
Q4. Your Scheme-C control series contains repeats and future cases. Explain, in terms of what the control series is supposed to represent, why deleting either would make the estimate worse, not better.
YOUR ANSWER:
6 Step 5 — One arithmetic, three estimands, in expectation (~8 min)
Your single draws landed at 2.46, 2.24, and 2.03. How much of that is the scheme, and how much is the luck of one draw? The chunk below (provided — just run it) repeats each scheme 500 times, then computes each scheme’s odds ratio from the pooled cells of all 500 replicates — the empirical version of the expected-cell-count argument from lecture.
set.seed(123)
R <- 500
rs_list <- lapply(seq_len(n_cases), function(j) {
which(coh$followup_days >= cases$followup_days[j] &
coh$patientunitstayid != cases$patientunitstayid[j])
})
one_rep <- function(scheme) {
if (scheme == "A") idx <- sample(nrow(pool_A), n_controls)
if (scheme == "B") idx <- sample(nrow(pool_B), n_controls)
if (scheme == "C") idx <- unlist(lapply(rs_list, function(rs)
rs[sample.int(length(rs), min(4, length(rs)))]))
pool <- switch(scheme, A = pool_A, B = pool_B, C = coh)
x <- sum(pool$early_vasopressor[idx]); n <- length(idx)
c(x = x, n = n, or = caco_or(exp_cases, n_cases, x, n))
}
reps <- bind_rows(lapply(c("A", "B", "C"), function(s)
t(replicate(R, one_rep(s))) |> as_tibble() |> mutate(scheme = s)))
pooled <- reps |>
group_by(scheme) |>
summarise(pooled_or = caco_or(exp_cases, n_cases, sum(x), sum(n)),
iqr_lo = quantile(or, .25), iqr_hi = quantile(or, .75))
pooled$target <- c(A = or_cohort, B = rr_cohort, C = ratr_cohort)[pooled$scheme]
pooled |> mutate(across(where(is.numeric), \(v) round(v, 3)))Checkpoint (pooled over 500 replicates vs. the cohort truth):
| scheme | pooled OR | target | middle half of single draws |
|---|---|---|---|
| A cumulative | 2.46 | odds ratio 2.46 | 2.05–3.09 |
| B case-base | 2.13 | risk ratio 2.13 | 1.88–2.46 |
| C incidence-density | 2.06 | rate ratio 2.00 | 1.74–2.73 |
targets <- tibble(scheme = c("A", "B", "C"),
target = c(or_cohort, rr_cohort, ratr_cohort))
reps |>
filter(is.finite(or)) |>
ggplot(aes(x = or)) +
geom_histogram(bins = 35, fill = "#56B4E9", color = "white") +
geom_vline(data = targets, aes(xintercept = target),
linetype = 2, linewidth = .6) +
facet_wrap(~ scheme, nrow = 1,
labeller = as_labeller(c(A = "A cumulative → cohort OR",
B = "B case-base → cohort RR",
C = "C density → cohort rate ratio"))) +
labs(x = "Exposure odds ratio from one case-control draw", y = "Count") +
theme_classic()This table and figure shows two things: First, schemes A and B sit exactly on their own targets, and each target is a different number — the estimand is set by the control-sampling scheme, not by the arithmetic. (Scheme C lands close to, though a shade above, the rate ratio: sampling a fixed 4 controls per case weights the late, tiny risk sets a little more than their person-time share; the fully time-matched analysis of Week 11 removes even that.) Second, look at the spread: any single case-control study — including yours from Steps 2–4 — is one draw from a wide distribution. Design fixes the target; sample size governs the noise.
Q5. Complete this sentence in your own words, then defend it in 3–4 more: “A case-control odds ratio is an estimator, not an estimand; what it estimates is decided by ______.”
YOUR ANSWER:
7 Step 6 — Save the samples for Week 11 (~2 min)
Week 11 returns to these exact datasets with regression tools (weighted logistic fits; conditional logistic for the time-matched series). Run the provided chunk to write them.
dir.create("data", showWarnings = FALSE)
write_csv(bind_rows(cases |> mutate(role = "case"),
ctrl_A |> mutate(role = "control")),
"data/caco_cumulative.csv")
write_csv(bind_rows(cases |> mutate(role = "case"),
ctrl_B |> mutate(role = "control")),
"data/caco_casebase.csv")
write_csv(bind_rows(cases |> mutate(role = "case"),
ctrl_C |> mutate(role = "control")),
"data/caco_density.csv")
length(list.files("data", pattern = "^caco_.*csv$")) # expect 3Checkpoint: 3 files in labs/lab2/data/.
8 Synthesis (~8 min)
The whole lab in one line per scheme:
| Scheme | Controls drawn from… | The exposure OR estimates… |
|---|---|---|
| Cumulative (survivor) | non-cases, at end of follow-up | cohort odds ratio |
| Case-base | the full cohort, at baseline | cohort risk ratio |
| Incidence-density | risk sets at each event time | cohort rate ratio |
Q6. A clinician asks you for “the 30-day risk of in-hospital death for an untreated patient like mine.” Explain why none of your three case-control samples can answer this on its own, and name two different additions that would restore an answer.
YOUR ANSWER:
Q7. In this lab the exposure was already recorded for all 1,274 patients, so sampling saved us nothing. Describe the two real-world conditions from lecture under which you would run one of these designs instead of analyzing the full cohort — and what today’s full-cohort benchmarks let us do that a real investigator never can.
YOUR ANSWER:
The takeaway. You ran the same cross-product on three samples from one cohort and got three systematically different answers — each one correct, for its own estimand. “What does a case-control odds ratio estimate?” has no answer until the control-sampling scheme is stated; once it is stated, the answer is exact, non-rare outcome and all. That is why the design section of a case-control paper is not boilerplate: it is the sentence that tells you what the number is.
Where this goes next: the Week 3 lecture supplies the algebra behind each scheme and extends the sampling idea to case-cohort designs and general outcome-dependent sampling; Week 11 returns to the three files you saved today, with regression machinery (weighted and conditional logistic) replacing the hand-computed cross-products.
References
- Lash TL, Rothman KJ. Case-Control Studies. In: Modern Epidemiology, 4th ed. (Lash TL, VanderWeele TJ, Haneuse S, Rothman KJ, eds.) Wolters Kluwer; 2021: chap. 8.
- Cornfield J. A method of estimating comparative rates from clinical data. J Natl Cancer Inst. 1951;11:1269–1275.
- Miettinen O. Estimability and estimation in case-referent studies. Am J Epidemiol. 1976;103(2):226–235.
- Greenland S. Model-based estimation of relative risks and other epidemiologic measures in studies of common outcomes and in case-control studies. Am J Epidemiol. 2004;160(4):301–305.
- O’Brien KM, Lawrence KG, Keil AP. The Case for Case-Cohort. Epidemiology. 2022;33(3):354–361.
- Pollard T, et al. eICU Collaborative Research Database Demo (v2.0.1). PhysioNet, 2021. doi:10.13026/4mxk-na84 (the Lab 1 source data).