Lab 2 — The Case-Cohort Design

EPID 785R · Regression and Study Design

Author

YOUR NAME HERE

Published

September 19, 2026

Introduction

Lab 1 focused on building an analytic dataset using target trial emulation concepts. From that lab, you should have a dataset with 1,292 ICU patients classified at a 60-minute landmark as early vasopressor initiators or non-initiators. Patients were followed from time to first ICU admission + 60 minutes (i.e., the landmark) to hospital discharge, with in-hospital death as the outcome. Covariates compiled in the analytic dataset include age, sex, baseline MAP, baseline heart rate, and the source from which baseline MAP was collected.

The file analytic_dataset.csv in the data folder represents this cohort constructed for a target trial emulation. We will use this file to construct a case-cohort dataset.

As a reminder (see Lecture 3 notes), we might carry out a case-cohort study if one of the covariates needed for the study was expensive to measure, if we had a dataset so large that our analysis is computationally inefficient, or if the outcome is rare.

The case-cohort design proceeds by extracting all the cases, and taking a random sample of everyone at baseline to create a subcohort. We then weight the sample back up using specific types of weights so that our analyses can be interpreted as being done in the full cohort.

In this lab, you will:

start by computing risks functions and contrasts in the original analytic dataset.
draw the case-cohort sample: a 20% subcohort plus all 136 cases.
construct an analytic case-cohort dataset in “long” (i.e., counting-process) form, and building the necessary weights for each row in this dataset. These weights should be a weight of “one” for each row representing the event occurrence, and a weight of 1/0.2 = 5 for all non-event rows.
we will then apply a weighted Kaplan-Meier estimator to compute risks, rates, the risk function, and a weighted Cox model for the hazard ratio; then a weighted log-Poisson regression to compute the risk ratio, and the weighted linear risk model for the risk difference. We’ll use the robust variance estimator to construct standard errors for these results.
Lastly, we’ll see what happens when we break the weights on purpose by keeping all the cases’ person-time (instead of only the random sub-cohort person time), and see what goes wrong when we do so.

Deliverables (to keep for yourself, since nothing is to be submitted): this document rendered to HTML with your answers written in the answer blocks, plus the case-cohort sample your code writes to data/ in Step 8. 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.Rproj in RStudio.
This lab needs your Lab 1 analytic dataset. But it is available to you in the data folder for lab2, included in the starter package on the course website.
R packages should include the tidyverse plus survival (Kaplan-Meier, Cox), and sandwich + lmtest (robust standard errors). You can use pacman::p_load(tidyverse, survival, sandwich, lmtest) to install anything missing.

pacman::p_load(tidyverse, survival, sandwich, lmtest)

## note: depending on where you put things on your computer, you may need to modify
## this code.

paths <- c("../lab1/data/derived/analytic_dataset.csv",  # Lab 1, completed in place
           "data/analytic_dataset.csv")                  # or the website 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 website into labs/lab2/data/.")
}
lab1 <- read_csv(path, show_col_types = FALSE)
nrow(lab1)  # expect 1,292

Instructions:

Chunks marked TODO contain ??? placeholders, which is invalid code, so the chunk will error out 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.

2 Step 1 — Estimates computed in the cohort directly (~8 min)

Before sampling anything from analytic_dataset.csv, let’s use it to establish the results 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). We don’t have time in this lecture/course to get into missing data methods, so we’ll conduct a complete case analysis, noting that this is not the ideal way to proceed with missingness here.

Complete TODO (1a) — the analytic cohort 2×2 table.

# 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 ICU stays, the 2×2 table is 6 / 21 (exposed) and 130 / 1,117 (unexposed). As a reminder from Lab 1: the unexposed group contains 107 patients who did actually start a vasopressor, but after the 60-minute window.

Now the benchmarks. Complete TODO (1b) (risks, risk ratio, risk difference) and TODO (1c) (person-time and rates; followup_days is each person’s time under observation, so summing it within exposure groups gives person-days at risk).

# Note: there are many ways to compute risks and rates, but we keep things
# simple here

# TODO (1b): risks, risk ratio, risk difference
R1 <- a1 / ???                       # risk in the exposed
R0 <- a0 / ???                       # risk in the unexposed
rr_cohort <- ??? / ???               # cohort risk ratio
rd_cohort <- ??? - ???               # cohort risk difference
round(c(R1 = R1, R0 = R0, RR = rr_cohort, RD = rd_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
round(c(pt1 = pt1, pt0 = pt0, total_pt = pt1 + pt0,
        rate1_per100pd = 100 * rate1, rate0_per100pd = 100 * rate0,
        rate_ratio = rate1 / rate0), 3)

Checkpoints: R1 = 0.222, R0 = 0.104; risk ratio = 2.13; risk difference = 0.118. Person-days 168.3 (exposed) and 7,310.1 (unexposed), 7,478.4 in total; rates 3.56 vs. 1.78 per 100 person-days; rate ratio = 2.00. Keep the total person-time in view — it returns as a design check in Step 3. (In Step 5 we will reproduce the risk ratio and risk difference with regression, which is how you will estimate them from the weighted sample.)

3 Step 2 — Draw the case-cohort sample (~8 min)

A case-cohort sample is a union of two pieces: a random subcohort drawn at baseline, and all the cases, whenever they occur. We use a 20% subcohort ($q = 0.20$): each cohort member flips the same $q$-coin, so every non-case in the sample is there with known probability $q$, and every case with probability 1. The expensive exposure is abstracted only for the union.

Complete TODO (2a): a chart gets abstracted if the person is in the subcohort or they became a case.

# There are many ways to draw a random sample from baseline. Here's one:

q <- 0.20                                  # subcohort sampling fraction
set.seed(13)
coh <- coh |>
  mutate(subc = rbinom(n(), 1, q) == 1)    # drawn using the rbinom function

# TODO (2a): the abstracted sample = subcohort OR case
samp <- coh |> filter(??? | ???)

c(n_subcohort   = sum(coh$subc),
  cases_inside  = sum(samp$hospital_death == 1 & samp$subc),
  cases_outside = sum(samp$hospital_death == 1 & !samp$subc),
  n_sample      = nrow(samp))

Checkpoints: 260 subcohort members; 136 cases, of whom 29 fell inside the subcohort by chance; 367 distinct people abstracted.

Note that if you changed the seed above, these numbers will not align with what you have.

The samp data frame now holds three types of people, and the whole analysis turns on handling them differently:

Type 1) non-events who were randomly sampled from the baseline data (subc = TRUE and hospital_death = 0);
Type 2) events who were randomly sampled from the baseline data (subc = TRUE and hospital_death = 1);
Type 3) events who were NOT sampled from the baseline data (subc = FALSE and hospital_death = 1).

Questions — write your answers below each prompt

Q1. Why must the subcohort be drawn blind to the outcome, and why is the 29-person overlap between the subcohort and the case series expected?

YOUR ANSWER:

Q2. Why can’t we simply analyze the 367 sampled people as a small cohort? That is, what, concretely, would be wrong with an estimate of the risk of death computed directly on this sample, without weights?

YOUR ANSWER:

4 Step 3 — Splitting up the person-time: the weights (~12 min)

Each of the three types of people in samp needs to be handled differently when constructing the weights:

Type	Person-time	Weight
Type 1 (subcohort non-event)	all their follow-up	$1/q = 5$
Type 2 (subcohort event)	split: start to right before the event; then the event	$1/q$, then $1$
Type 3 (event outside the subcohort)	split: start to right before the event gets weight 0 (that person-time is removed); then the event	—, then $1$

Why do the weights vary over time? Because a case’s probability of being in the sample changes at their event: before it, they could only enter through the subcohort (probability $q$); at the event, they are sampled with probability 1. Splitting the person-time lets each segment carry the right inverse-probability weight. Note carefully: this is not about time-varying exposures or covariates — the exposure is fixed at baseline. The split exists because the weight changes over time. And a Type-3 case contributes no comparison person-time at all: the subcohort, scaled by $1/q$, already represents everyone’s person-time — including the follow-up of people who will eventually become cases — so their earlier follow-up is not used, because it is already counted.

We encode the split in counting-process (“long”) format: one row per person-time segment, (start, stop], with an event indicator and a weight w. eps separates “right before the event” from the event itself; our event times sit on a 1-minute grid (1/1440 of a day), so eps = 1e-4 days (about 9 seconds) slips below it.

Complete TODO (3a) (where a Type-2 comparison row stops) and TODO (3b) (the event rows: where they start, and their weight).

eps <- 1e-4   # smaller than the 1-minute grid the event times sit on

cc <- rbind(
  # Type 1 — subcohort non-events: full follow-up at weight 1/q
  with(subset(samp, subc & hospital_death == 0),
       data.frame(id = patientunitstayid, a = early_vasopressor,
                  start = 0, stop = followup_days,
                  event = 0, w = 1/q)),
  # Type 2, pre-event segment — subcohort events: weight 1/q from start
  # to right before the event
  with(subset(samp, subc & hospital_death == 1),
       data.frame(id = patientunitstayid, a = early_vasopressor,
                  # TODO (3a): this comparison row stops just BEFORE the event
                  start = 0, stop = ??? - eps,
                  event = 0, w = 1/q)),
  # Types 2 & 3, event segment — every case's event at weight 1
  # (Type 3's pre-event person-time gets weight 0: it has no row at all)
  with(subset(samp, hospital_death == 1),
       data.frame(id = patientunitstayid, a = early_vasopressor,
                  # TODO (3b): the event row covers only the sliver just
                  # before the event, at weight 1
                  start = pmax(0, ??? - eps),
                  stop = followup_days, event = 1, w = ???))
)

# rows per type: Type 1 = 231; Type 2 = 29 people x 2 rows = 58; Type 3 = 107
nrow(cc)   # expect 231 + 58 + 107 = 396 rows

Checkpoint: 396 rows from 367 people — each of the 29 Type-2 people produces two rows; each Type-1 and Type-3 person produces one.

Design checks. Before estimating anything, verify that the weighted split reconstructs the cohort — step 3 of the lecture’s implementation checklist. Weighted person-time should approximate the cohort’s total (it is an estimate built from a 20% sample); the weighted event count should match the cohort’s 136 exactly (every case is in, at weight 1).

c(cohort_pt   = sum(coh$followup_days),
  weighted_pt = sum(cc$w * (cc$stop - cc$start)))   # expect ~7,478 vs ~7,665
sum(cc$w[cc$event == 1])    # weighted events: expect exactly 136
sum(cc$w[cc$start == 0])    # weighted size of the time-zero risk set: expect 1,300

Checkpoints: weighted person-time 7,664.6 against the cohort’s 7,478.4 (off by 2.5% — a sampling estimate doing its job); weighted events exactly 136; the weighted risk set at time zero is 1,300 ≈ 1,274 people.

Questions

Q3. Why must a Type-2 person’s comparison row stop just before their event rather than at it? Say precisely what would be counted twice if it ran to the event at weight 5 — and connect it to lecture’s warning about the subcohort’s cases.

YOUR ANSWER:

5 Step 4 — Time-to-event analyses: weighted KM and Cox (~15 min)

Kaplan-Meier and Cox are built on risk sets indexed by time, so they need the person-time split: the weight a person carries must be allowed to change over their follow-up — $1/q$ while serving as comparison person-time, 1 at the event, 0 (no row) for the never-sampled comparison time of Type-3 cases. That is exactly what cc encodes.

First, the risk function. The Kaplan-Meier estimator turns risk sets and events into a cumulative risk curve $F(t) = 1 - S(t)$; fed the weighted rows, each subcohort member holds a place in every risk set for five people, and each case’s event arrives once, at weight 1. We fit it in the case-cohort sample, then run the same analysis in the original cohort, and overlay the two. Run the provided chunk.

km_cc   <- survfit(Surv(start, stop, event) ~ 1, data = cc, weights = w)
km_full <- survfit(Surv(followup_days, hospital_death) ~ 1, data = coh)

km_df <- bind_rows(
  tibble(time = km_full$time, risk = 1 - km_full$surv,
         analysis = "Full cohort (1,274 charts)"),
  tibble(time = km_cc$time, risk = 1 - km_cc$surv,
         analysis = "Case-cohort, weighted (367 charts)"))

km_df |>
  filter(time <= 30) |>
  ggplot(aes(x = time, y = risk, color = analysis)) +
  geom_step(linewidth = .7) +
  scale_color_manual(values = c("Full cohort (1,274 charts)" = "grey25",
                                "Case-cohort, weighted (367 charts)" = "#56B4E9")) +
  labs(x = "Days since the landmark", y = "Cumulative risk of in-hospital death",
       color = NULL) +
  theme_classic() + theme(legend.position = "bottom")

# read both curves at a few horizons
risk_at <- function(fit, t) 1 - summary(fit, times = t)$surv
round(rbind(full_cohort = risk_at(km_full, c(5, 10, 30)),
            case_cohort = risk_at(km_cc, c(5, 10, 30))), 3)

Checkpoints: the two step functions ride on top of each other; at 5, 10, and 30 days the full-cohort risks are 0.099 / 0.154 / 0.335 and the weighted case-cohort risks 0.096 / 0.151 / 0.354. Ten-day risk of in-hospital death: 15.4% from 1,274 charts, 15.1% from 367.

One honest footnote on what these curves are: Kaplan-Meier here treats live discharge as censoring — as if discharged patients remained at risk of in-hospital death — which is why the 30-day figure (0.335) towers over the crude 136/1,274 = 0.107. Whether that is a sensible estimand is a real question, and it is Week 12’s question (censoring, competing events). Today the curves serve a different purpose: the same estimator, fed the weighted sample, reproduces whatever the full cohort would have said. The design is validated on its own terms.

Next, the weighted Cox model — the estimator lecture developed the case-cohort design around, and cc is already in exactly the right shape. We fit it twice in the case-cohort sample — once accepting the default model-based (naive) standard error, once with the robust variance clustered on the person, so a Type-2 person’s two rows count as one individual — and then fit the same model in the original cohort.

cox_cc_naive  <- coxph(Surv(start, stop, event) ~ a,
                       data = cc, weights = w, robust = FALSE)
cox_cc_robust <- coxph(Surv(start, stop, event) ~ a,
                       data = cc, weights = w, cluster = id)
cox_full <- coxph(Surv(followup_days, hospital_death) ~ early_vasopressor,
                  data = coh)

round(rbind(
  full_cohort  = c(loghr = coef(cox_full), se = sqrt(vcov(cox_full)[1, 1])),
  cc_naive_SE  = c(coef(cox_cc_naive),  sqrt(cox_cc_naive$var[1, 1])),
  cc_robust_SE = c(coef(cox_cc_robust), sqrt(cox_cc_robust$var[1, 1]))), 3)

Checkpoints: full-cohort log hazard ratio 0.722 (SE 0.419); the case-cohort fit lands at 0.725 — on top of the full cohort — with naive SE 0.419 and robust SE 0.655.

Look hard at that naive SE: it is the full cohort’s standard error, to three decimals, and coxph() hands it to you without complaint. That is not luck: the weighted risk sets are scaled back up to full-cohort size, so a model-based variance formula “sees” 1,274 people and reports the precision that 1,274 abstracted charts would have earned. We abstracted 367. The robust, person-clustered standard error charges for the difference — had a different 20% subcohort been drawn, the risk sets (and the estimate) would have differed. Report the robust one, always.

Questions

Q4. Explain, mechanically, why the naive SE of the weighted Cox fit must reproduce the full-cohort SE — and why that precision is a claim this study design cannot back up. Where, physically, did the missing information go?

YOUR ANSWER:

Q5. Write the sentence you would put in a methods section reporting the robust Cox analysis: point estimate, interval, variance estimator, and the design details a reader needs to reproduce it.

YOUR ANSWER:

Optional — the curves by exposure arm

Refit both survfit() calls with ~ early_vasopressor (full cohort) and ~ a (case-cohort) and overlay the four curves. Two things to notice, one per design lesson. First, the weighted exposed-arm curve is carried by about 10 people — 6 cases plus 4 subcohort non-cases — so it is steppy and fragile: with 27 exposed patients in the whole cohort, no 20% subcohort supports a smooth exposed-arm curve (Step 7 quantifies this). Second, in the full cohort the exposed curve leaps early (0.216 by day 5), then runs out of people — the last exposed patient leaves at day 16, freezing the curve at 0.276 — while the unexposed curve climbs past it (0.312 by day 20, 0.334 by day 30). Exposed deaths concentrated early plus discharge treated as censoring: curve-crossing here is exactly the kind of artifact Week 12 teaches you to interrogate rather than interpret away.

6 Step 5 — Binary-outcome analyses: one weight per person (~12 min)

The log-Poisson (risk ratio) and linear (risk difference) models are fit to one binary outcome per person. Person-time never enters these models explicitly, so there is nothing to split: each person can carry a single weight — the weight of their final person-time segment from Step 3. Type 1 → $1/q$ (their whole follow-up is one segment at $1/q$); Types 2 and 3 → 1 (their final segment is the event, at weight 1). Same weighting logic as Step 4, collapsed to one number per person because these models have no time axis.

Complete TODO (5a): extract each person’s final segment from cc.

# TODO (5a): each person's final segment = the row with their largest
# value of which variable?
w_exit <- cc |>
  group_by(id) |>
  slice_max(???, n = 1, with_ties = FALSE) |>
  ungroup() |>
  select(patientunitstayid = id, w)

samp <- samp |> left_join(w_exit, by = "patientunitstayid")

samp |> select(patientunitstayid, followup_days, death_event, w, subc)

# the weighted 2x2, next to the cohort's
samp |>
  count(early_vasopressor, hospital_death, wt = w, name = "weighted_n") |>
  mutate(cohort_n = c(b0, a0, b1, a1))
sum(samp$w)   # the weighted sample "size"

Checkpoints: scrolling the printed rows, every event carries w = 1 and every sampled non-event w = 5. The weighted 2×2 is 6 / 20 (exposed) and 130 / 1,135 (unexposed), sitting next to the cohort’s 6 / 21 and 130 / 1,117; the weighted sample size is 1,291 ≈ 1,274. The case cells are recovered exactly (every case is in, at weight 1); the non-case cells are estimates with sampling variability (4 sampled exposed non-cases × 5 = 20, standing in for the true 21).

Now the two effect measures, by weighted regression. A log-link Poisson model of a binary outcome (the “modified Poisson” procedure from Zou 2004 AJE) exponentiates to the risk ratio; an identity-link linear model of the same outcome returns the risk difference. The chunk fits each in the case-cohort sample (with weights = w), then fits the same two models in the original cohort.

# case-cohort, weighted
fit_rr_cc <- glm(hospital_death ~ early_vasopressor,
                 family = poisson(link = "log"),
                 weights = w, data = samp)
summary(fit_rr_cc)$coefficients                       # naive (model-based) SE
coeftest(fit_rr_cc, vcov = vcovHC(fit_rr_cc, type = "HC3"))  # robust SE
exp(coef(fit_rr_cc)["early_vasopressor"])             # the estimated risk ratio

fit_rd_cc <- lm(hospital_death ~ early_vasopressor,
                weights = w, data = samp)
coeftest(fit_rd_cc, vcov = vcovHC(fit_rd_cc, type = "HC3"))  # robust only — see below

# the same models in the original cohort
fit_rr_full <- glm(hospital_death ~ early_vasopressor,
                   family = poisson(link = "log"), data = coh)
coeftest(fit_rr_full, vcov = vcovHC(fit_rr_full, type = "HC3"))
exp(coef(fit_rr_full)["early_vasopressor"])   # the cohort risk ratio, again

fit_rd_full <- lm(hospital_death ~ early_vasopressor, data = coh)
coeftest(fit_rd_full, vcov = vcovHC(fit_rd_full, type = "HC3"))

And side by side:

glm_row <- function(fit, exponentiate = FALSE) {
  b  <- unname(coef(fit)["early_vasopressor"])
  nv <- unname(summary(fit)$coefficients["early_vasopressor", "Std. Error"])
  rb <- unname(sqrt(vcovHC(fit, type = "HC3")["early_vasopressor",
                                              "early_vasopressor"]))
  est <- c(b, b - 1.96 * rb, b + 1.96 * rb)
  if (exponentiate) est <- exp(est)
  c(estimate = est[1], naive_se = nv, robust_se = rb,
    lcl = est[2], ucl = est[3])
}

# risk ratio (estimate and CI on the RR scale; SEs on the log scale)
round(rbind(full_cohort = glm_row(fit_rr_full, exponentiate = TRUE),
            case_cohort = glm_row(fit_rr_cc,   exponentiate = TRUE)), 3)

# risk difference (robust CI; naive column dropped — robust only for the RD)
round(rbind(full_cohort = glm_row(fit_rd_full),
            case_cohort = glm_row(fit_rd_cc))[, -2], 3)

Checkpoints:

	estimate	naive SE	robust SE	95% CI
risk ratio, full cohort	2.132	0.418	0.383	1.006–4.516
risk ratio, case-cohort	2.246	0.418	0.586	0.712–7.079
risk difference, full cohort	0.118	—	0.084	−0.046–0.282
risk difference, case-cohort	0.128	—	0.134	−0.134–0.390

Notice something you have already seen once today, in Step 4’s Cox table: the two naive SEs in the risk-ratio table are identical (0.418 and 0.418) even though one analysis measured 1,274 exposures and the other 367 — the weighted pseudo-population has size $\sum_i w_i \approx 1{,}274$, so the model-based variance formula reports full-cohort precision the abstraction budget never bought. The robust column tells the truth, widening from 0.383 to 0.586. (For the risk difference we report only the robust SE; a weighted linear model’s naive SE partially absorbs the sampling inflation through its residual variance — 0.114 here, much of the way to robust — so the naive-vs-robust contrast is muddier than the Poisson’s. The rule that is clean: sampling weights → robust SE, always.)

Questions

Q6. In the weighted 2×2, which cells are exact and which are estimates? Connect this to the two pieces of the design.

YOUR ANSWER:

7 Step 6 — Break it on purpose (~10 min)

This step is an illustrative example, separate from the main lesson: do not do this in a real analysis.

Here is the mistake this design invites: “We paid to abstract all 136 cases — why throw their follow-up away? Keep every case’s full person-time.” The chunk below rebuilds the analytic rows with exactly one change: each case’s row now runs from 0 to their event (start = 0), at weight 1, instead of covering only the final sliver. The subcohort comparison rows are untouched. Run it and re-run the design checks, the KM readouts, and the Cox model.

cc_bad <- rbind(
  cc[cc$event == 0, ],    # subcohort comparison rows, unchanged
  with(subset(samp, hospital_death == 1),
       data.frame(id = patientunitstayid, a = early_vasopressor,
                  start = 0,                     # <- the only change
                  stop = followup_days, event = 1, w = 1))
)

# the design check catches it: person-time inflated, events unchanged
c(weighted_pt_bad = sum(cc_bad$w * (cc_bad$stop - cc_bad$start)),
  weighted_events_bad = sum(cc_bad$w[cc_bad$event == 1]))

# effect on the KM: risks depressed at every horizon
km_bad <- survfit(Surv(start, stop, event) ~ 1, data = cc_bad, weights = w)
round(rbind(full_cohort    = risk_at(km_full, c(5, 10, 30)),
            case_cohort    = risk_at(km_cc, c(5, 10, 30)),
            double_counted = risk_at(km_bad, c(5, 10, 30))), 3)

# effect on the overall death rate (per 100 person-days)
round(100 * c(
  full_cohort = (a1 + a0) / sum(coh$followup_days),
  case_cohort = sum(cc$w[cc$event == 1]) / sum(cc$w * (cc$stop - cc$start)),
  double_counted = sum(cc_bad$w[cc_bad$event == 1]) /
    sum(cc_bad$w * (cc_bad$stop - cc_bad$start))
), 3)

# effect on the Cox model
cox_bad <- coxph(Surv(start, stop, event) ~ a,
                 data = cc_bad, weights = w, cluster = id)
round(rbind(
  full_cohort    = c(loghr = coef(cox_full), se = sqrt(vcov(cox_full)[1, 1])),
  case_cohort    = c(coef(cox_cc_robust), sqrt(cox_cc_robust$var[1, 1])),
  double_counted = c(coef(cox_bad), sqrt(cox_bad$var[1, 1]))), 3)

Checkpoints: the event check still passes (136), but the person-time check now fails loudly: 8,397.1 person-days against the cohort’s 7,478.4 — inflated by 732.5, which is exactly the 136 cases’ total follow-up, now counted on top of the subcohort’s $1/q$-weighted representation of it. Every risk the double-counted curve reports is too low (0.089 / 0.139 / 0.324 at days 5 / 10 / 30); the overall death rate reads 1.62 per 100 person-days against the correct 1.77 and the cohort’s 1.82 — about 9% low, at every horizon, in every seed. The Cox log hazard ratio, though, barely moves: 0.716 against the correct 0.725.

Add the double-counted curve to the overlay and look at it:

km_df3 <- bind_rows(
  km_df,
  tibble(time = km_bad$time, risk = 1 - km_bad$surv,
         analysis = "Double-counted: all case person-time kept"))

km_df3 |>
  filter(time <= 30) |>
  ggplot(aes(x = time, y = risk, color = analysis)) +
  geom_step(linewidth = .7) +
  scale_color_manual(values = c(
    "Full cohort (1,274 charts)" = "grey25",
    "Case-cohort, weighted (367 charts)" = "#56B4E9",
    "Double-counted: all case person-time kept" = "#D55E00")) +
  labs(x = "Days since the landmark", y = "Cumulative risk of in-hospital death",
       color = NULL) +
  theme_classic() + theme(legend.position = "bottom")

The double-counted curve is not noisy — it is systematically depressed, sitting at about 92% of the correct curve everywhere. The mechanism: cases are in the sample because they are cases, so their person-time is an outcome-selected sample of the cohort’s person-time. No weight attached to it can fix that; the only valid move is the one lecture prescribes — their pre-event follow-up “will not be used as comparison person-time, because the subcohort already represents it.”

One more reading of the damage. The ratio measures barely notice: both arms’ denominators inflate by similar fractions (8.3% exposed, 9.8% unexposed), so the double-counted rate ratio and hazard ratio land near 2.0 anyway. If you only ever report ratios, this error can hide indefinitely. But absolute risks and rates — the very quantities the case-cohort design, unlike case-control, promises to deliver — are wrong by ~10%, everywhere.

An aside — the same mistake where the truth is known

Is the unmoved hazard ratio a property of the mistake, or luck? We ran the identical error in a simulated world where the answer is known by construction: a cohort of 100,000 with 50% exposure, exponential event times with a true hazard ratio of 2, ~52% of the cohort experiencing the event by the end of follow-up, and the same $q = 0.20$ case-cohort construction as this lab. The result:

analysis	HR
truth	2.00
full cohort	1.97
case-cohort, correct split	2.04
case-cohort, double-counted	1.79

With a common outcome there is far more case person-time to double count (the person-time check overshoots by +33% there, versus +12% here), and the exposed arm — where events concentrate — retains proportionally more of it: the same one-token error that barely moved our hazard ratio drags the simulated one about 10 standard errors below the truth, toward the null. The ratio’s survival in our data was luck; the person-time design check catches the error in both worlds.

Questions

Q7. The event check passed while the person-time check failed, by exactly 732.5 person-days. Whose time was double-counted, and by which two routes does it enter the double-counted denominator?

YOUR ANSWER:

Q8. Your colleague replies: “But the double-counted rate ratio and hazard ratio were basically right — 2.0! So the shortcut is fine.” Give the two-sentence rebuttal: why the ratios survived in this cohort, and why that is not exoneration.

YOUR ANSWER:

8 Step 7 — One draw vs. the design, in expectation (~7 min)

Your Step-5 estimates (RR 2.25, RD 0.128) sit near the truth (2.13, 0.118) but not on it. How much of the gap is the design, and how much is the luck of one 20% draw? The chunk below (provided — just run it) redraws the subcohort 2,000 times, computing each draw’s weighted risk ratio, risk difference, and death rate from the weighted cells directly, plus the double-counted rate from Step 6’s broken person-time. Pooling across draws gives the design’s expected behavior, next to the full-cohort truth.

set.seed(42)
R    <- 2000
fu   <- coh$followup_days
expo <- coh$early_vasopressor
dth  <- coh$hospital_death
case_pt <- sum(fu[dth == 1])

one_rep <- function() {
  subc <- rbinom(length(fu), 1, q) == 1
  k1 <- sum(subc & expo == 1 & dth == 0)   # exposed non-cases sampled
  k0 <- sum(subc & expo == 0 & dth == 0)   # unexposed non-cases sampled
  # weighted person-time: subcohort comparison rows, then the case rows
  pt_sub <- (sum(fu[subc & dth == 0]) +
               sum(pmax(fu[subc & dth == 1] - eps, 0))) / q
  slivers <- sum(fu[dth == 1] - pmax(fu[dth == 1] - eps, 0))
  c(k1 = k1, k0 = k0,
    rr = (a1 / (a1 + k1 / q)) / (a0 / (a0 + k0 / q)),
    rd = a1 / (a1 + k1 / q) - a0 / (a0 + k0 / q),
    pt_w   = pt_sub + slivers,   # correct: cases enter as slivers
    pt_bad = pt_sub + case_pt)   # double-counted: full case follow-up
}

reps <- as_tibble(t(replicate(R, one_rep())))

## Pool the weighted cells / person-time across draws and set each
## quantity next to its full-cohort value. The iqr columns show the
## middle half of the SINGLE-draw estimates: wide for the ratio (noise),
## while the double-counted rate is narrow but off target (bias).
pooled_rr   <- (a1 / (a1 + mean(reps$k1) / q)) / (a0 / (a0 + mean(reps$k0) / q))
pooled_rd   <- a1 / (a1 + mean(reps$k1) / q) - a0 / (a0 + mean(reps$k0) / q)
cohort_rate <- 100 * (a1 + a0) / (pt1 + pt0)
rate_w      <- 100 * 136 / reps$pt_w     # each draw's correctly weighted rate
rate_bad    <- 100 * 136 / reps$pt_bad   # each draw's double-counted rate

iqr <- function(x) c(iqr_lo = unname(quantile(x, .25)),
                     iqr_hi = unname(quantile(x, .75)))

round(rbind(
  risk_ratio          = c(full_cohort = rr_cohort,  pooled_draws = pooled_rr,
                          iqr(reps$rr)),
  risk_difference     = c(rd_cohort,   pooled_rd,   iqr(reps$rd)),
  rate_per100pd       = c(cohort_rate, 100 * 136 / mean(reps$pt_w),   iqr(rate_w)),
  rate_double_counted = c(cohort_rate, 100 * 136 / mean(reps$pt_bad), iqr(rate_bad))
), 3)

Checkpoint (2,000 subcohort draws vs. the full cohort):

quantity	full cohort	pooled over draws	middle half of single draws
risk ratio	2.132	2.125	1.727–2.791
risk difference	0.118	0.117	0.081–0.183
death rate (per 100 pd)	1.819	1.819	1.720–1.935
double-counted rate	1.819	1.657	1.574–1.752

Two lessons in one table. The design is valid: in expectation, the weighted sample recovers the cohort’s risk ratio, risk difference, and rate essentially exactly — while the double-counted rate misses its target in every draw (its entire middle half sits below 1.819): bias, not noise. And the design is honest about its price: any single 20% draw of this cohort estimates the risk ratio with real spread, because the exposed non-case cell rests on a handful of people (about 4, in expectation — your seed-13 draw got exactly 4).

Questions

Q9. The single-draw spread is driven almost entirely by one random quantity. Name it, explain why it is so influential in this cohort, and explain why the obvious fix — oversample exposed people into the subcohort — cannot be implemented as stated. What could an investigator legitimately dial instead?

YOUR ANSWER:

9 Step 8 — Save the sample for Week 11 (~2 min)

Week 11 returns to outcome-dependent designs with the full regression toolkit, and reuses this sample. Run the provided chunk.

dir.create("data", showWarnings = FALSE)
write_csv(samp, "data/case_cohort_sample.csv")
nrow(samp)   # expect 367 rows, with subc and w attached

Checkpoint: data/case_cohort_sample.csv, 367 rows, 15 columns (the 13 Lab 1 columns plus subc and w).

10 Optional challenge — turn the design dial (~5 min)

Optional

The subcohort fraction $q$ is the design’s cost–precision dial. Change q in Step 2, re-run Steps 2, 3, and 5 (the seed keeps each draw reproducible per q), and record the robust SE of the log risk ratio; then restore q <- 0.20 and re-run before rendering. You should see, approximately:

$q$	charts abstracted	cost	naive SE	robust SE
0.10	261	$39,150	0.418	0.949
0.20	367	$55,050	0.418	0.586
0.40	584	$87,600	0.418	0.448

(The full cohort abstracts 1,274 charts at $191,100 for robust SE 0.383.) The point estimates wobble from draw to draw — that is Step 7’s lesson — but the standard-error columns are the design speaking: the naive column never moves, no matter how little you abstract, which is everything wrong with it in one column; the robust column pays for exactly what you didn’t measure, and buys precision back as $q$ rises.

11 Synthesis (~7 min)

You have now executed, on real data, the applied checklist from lecture’s case-cohort section:

Lecture checklist	Where you did it
1. Construct the sample (subcohort of fraction $q$ + all cases)	Step 2
2. Define inclusion indicators and weights (time-varying; $1/q$ then 1)	Step 3 (split), Step 5 (collapsed to one per person)
3. Check the weighted sample against the cohort	Step 3 (and watched it fail in Step 6)
4. Fit the weighted regressions (Cox; GLMs for risk models)	Steps 4–5
5. Compute the robust variance, clustered on the person	Steps 4–5
6. Report the estimates with the design	Q5

The takeaway. A 20% subcohort plus the cases — 367 charts out of 1,274 — recovered the cohort’s risk curve, hazard ratio, risk ratio, and risk difference, because every analysis step “knew” the sampling design through the weights: the time-to-event analyses through the split person-time ($1/q$ for the subcohort’s segments, 1 for each event, and nothing at all for case person-time the subcohort already represents), and the binary-outcome regressions through each person’s final-segment weight. The two ways this analysis dies quietly are now both in your hands: a model-based standard error that reports precision the abstraction budget never bought, and outcome-selected person-time smuggled into the denominators — caught, respectively, by the robust variance and by the design check that weighted person-time must match the cohort.

Where this goes next: Week 4 begins the regression toolkit that Steps 1 and 5 borrowed from (glm, links, and what “saturated” means); Week 7 opens up the sandwich variance; Week 11 returns to outcome-dependent sampling with the full toolkit — including the case-cohort sample you saved in Step 8 — and Weeks 12–13 give the risk curves and hazard ratios of Step 4 their honest treatment of censoring and competing events.

References

O’Brien KM, Lawrence KG, Keil AP. The Case for Case-Cohort: An Applied Epidemiologist’s Guide to Reframing Case-Cohort Studies to Improve Usability and Flexibility. Epidemiology. 2022;33(3):354–361.
Prentice RL. A case-cohort design for epidemiologic cohort studies and disease prevention trials. Biometrika. 1986;73:1–11.
Barlow WE. Robust variance estimation for the case-cohort design. Biometrics. 1994;50:1064–1072.
Therneau TM, Li H. Computing the Cox model for case cohort designs. Lifetime Data Anal. 1999;5:99–112.
Zou G. A modified Poisson regression approach to prospective studies with binary data. Am J Epidemiol. 2004;159(7):702–706.
Pollard T, et al. eICU Collaborative Research Database Demo (v2.0.1). PhysioNet, 2021. doi:10.13026/4mxk-na84 (the Lab 1 source data).

Type	Person-time	Weight
Type 1 (subcohort non-event)	all their follow-up	\(1/q = 5\)
Type 2 (subcohort event)	split: start to right before the event; then the event	\(1/q\), then \(1\)
Type 3 (event outside the subcohort)	split: start to right before the event gets weight 0 (that person-time is removed); then the event	—, then \(1\)