Lab 3 — Three Questions

EPID 785R · Regression and Study Design

Author

YOUR NAME HERE

Published

September 27, 2026

Introduction

In lab 1, we 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. In lab 2, we treated that cohort as a fully enumerated source population and sampled from it to conduct a case-cohort analysis.

In this lab, we’ll take a bit of a different approach and analyze the ICU cohort using regression, and we’ll use this regression model for three different kinds of questions: a descriptive one, a predictive one, and a causal one.

The goal of this lab will be to develop an understanding of what kinds of considerations are entailed by descriptive analyses, predictive analyses, and causal analyses. Of course, we don’t have time to work through the details of each of these problem areas. Rather, what we’ll try to do is connect this practical analysis to some of the concepts that we covered in lecture, namely: when we use the machinery of regression, it is fundamentally agnostic to questions of description, prediction, or causal inference.

What changes across the three questions is everything outside the model, such as the target that we want to quantify (the estimand), and the assumptions needed to connect the fitted output to it.

In this lab, you will:

fit the model an analyst might fit when no question has been asked yet, and read its printout;
write out a descriptimand (following Lesko, Fox, and Edwards’ framework) and decide what the printout can and cannot say about it;
write a predictimand (following Vansteelandt and Steen) and evaluate the model as a forecasting tool;
revisit the causal estimand you formulated in Lab 1 and price the two layers of assumptions the printout would need to answer it;
step back and notice what changed across the three readings (the question, the estimand, the assumptions) and what did not (the printout).

Deliverables (yours to keep, nothing is submitted): this document rendered to HTML with your answers in the answer blocks. Only two chunks contain code TODOs; most of the blanks in this lab are sentences, because most of the work regression doesn’t do for you is done in sentences. Once the two TODO chunks run, flip eval: false to eval: true in the YAML header and re-render.

Background reading for this lab includes: Carlin & Moreno-Betancur (2025) and the Vansteelandt & Steen commentary (both in the same PDF on the course site), plus Lesko, Fox & Edwards (2022) for Step 1.

1 Before you start (~4 min)

Open lab3.Rproj in RStudio.
The lab needs your Lab 1 analytic dataset. As in Lab 2, a copy is also included in this lab’s data/ folder in the starter package.
Packages: just the tidyverse this time.

pacman::p_load(tidyverse)

## 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 the website into labs/lab3/data/.")
}
lab1 <- read_csv(path, show_col_types = FALSE)
nrow(lab1)   # expect 1,292

# As in Lab 2: 18 stays have a blank hospital discharge status. We again
# proceed with a complete-case analysis, again noting that this is for 
# convenience.
coh <- lab1 |> filter(!is.na(hospital_death))
nrow(coh)    # expect 1,274

2 Step 0 — A model before a question (~4 min)

Here is how an enormous amount of empirical research actually begins: with a dataset, an outcome, and the covariates that happen to be available, all handed to a fitting routine before anyone has said what the analysis is for.

Run the chunk (it is fully provided).

fit <- glm(hospital_death ~ early_vasopressor + age + sex +
             baseline_map + baseline_hr,
           family = binomial, data = coh)

round(summary(fit)$coefficients, 3)   # the printout
round(exp(coef(fit)), 3)              # the same, as odds ratios
nobs(fit)                             # how many rows did the fit use?

Checkpoints: the coefficient table reads, to three decimals:

term	estimate	SE	p	OR
(Intercept)	−6.877	1.083	<0.001	—
early_vasopressor	0.429	0.509	0.400	1.54
age (per year)	0.046	0.007	<0.001	1.047
sexMale	0.290	0.195	0.137	1.34
baseline_map (per mmHg)	−0.012	0.013	0.339	0.99
baseline_hr (per bpm)	0.024	0.004	<0.001	1.024

And nobs(fit) says 1,271 — though coh has 1,274 rows. Three people were removed in the fitting process: the patients had missing heart rate at baseline (Lab 1’s QC report counted them). glm(), by default, cannot handle any missing data, so they are dropped silently.

We’ll use printouts like this, as well as functions of printouts like this, as the main output for the lab. The next three steps ask it, in turn, a descriptive, a predictive, and a causal question.

We’ll see that the numbers don’t change (they can’t change), but their meaning, and the assumptions needed to sustain that meaning, change completely.

3 Step 1 — Reading it descriptively (~12 min)

Start with the question in its usual, ill-posed form:

“What’s the mortality of hypotensive ICU patients?”

Lesko, Fox, and Edwards argue that a well-defined descriptive question must state: (1) the target population, characterized by person and place and anchored in time; (2) the outcome or health state of specific interest; (3) the measure of occurrence used to summarize the outcome; and (4) any auxiliary variables, with their roles (stratification vs. standardization).

Write the specific, sharpened descriptive question for our setting. Everything you need is in the Lab 1 protocol: Who is in the cohort? What is the anchor? Where do the data come from (the eICU Collaborative Research Database Demo: US hospitals participating in a tele-ICU program, 2014–2015)? What is the outcome? What summary measure do we want to use?

Your descriptimand — complete the template

Among ________________________________ (person — who, anchored when?) in ________________________________ (place) during ________ (time), the ________________________ (measure of occurrence) of ________________________ (outcome).

YOUR ANSWER:

Now estimate it. Complete TODO (1a) — the “null model” from lecture: an intercept-only regression, which is the simplest descriptive machine there is.

# the estimate, computed as a proportion
sum(coh$hospital_death) / nrow(coh)

# TODO (1a): the same estimate, computed by regression — a model with an
# intercept and nothing else (in R's formula language, "~ 1")
fit0 <- glm(hospital_death ~ ???, family = binomial(link = "logit"), data = coh)
plogis(coef(fit0))    # inverse-logit of the intercept: you should know what this equation is
nobs(fit0)

# a prespecified stratification: the outcome distribution by sex
coh |>
  group_by(sex) |>
  summarise(n = n(), deaths = sum(hospital_death),
            risk = round(mean(hospital_death), 3))

Checkpoints: the proportion is 136/1,274 = 0.107, and plogis(coef(fit0)) returns exactly the same 0.107. The intercept-only model targets the sample proportion of the outcome, using the machinery of regression. Note that we could have done this with other non-regression tools as well. Its nobs is 1,274: this model did not remove the individuals with missing heart rate data, because heart rate is not in the model. By sex: 86/754 = 0.114 among males, 50/520 = 0.096 among females (a crude odds ratio of about 1.21, if you compute it).

In the print-out from Step 0, the coefficient for sexMale yielded an OR = 1.34. The crude comparison you just made says 1.21. Why are they different? Because they answer different questions. The OR = 1.21 value compares the men in this cohort to the women in this cohort. However, the OR = 1.34 compares men to women of the same age, baseline MAP, and heart rate. This is a statement about strata that our descriptive question (at least as written) did not specify. The takeaway here, as noted in Lesko et al. is that adjustment in a descriptive analysis “implies an intervention on the data”. In this case, it’s a reweighting toward a population that, in this case, was never named. If the descriptimand is the sex contrast in this population, our goal should be to target the crude number, not the adjusted estimate.

Another important consideration is that the estimate of 0.107 characterizes the analytical sample. Whether 0.107 describes the target population of interest depends on assumptions the model cannot see: for instance, that the eICU demo’s patients stand in for that population (hospitals chose to join a tele-ICU program, and the demo is a non-random subsample of the full database, which is itself a convenience sample of units using one vendor’s tele-ICU platform, not a random sample of any population of interest). Additionally, these data come from measurement devices which, themselves, may be subject to error. Sampling and measurement therefore become key to understanding the extent to which we can interpret any result from a regression model as a “descriptimand” of interest.

Questions — write your answers below each prompt

Q1. In Lab 2, a Kaplan-Meier analysis of this same cohort reported a 30-day “risk of in-hospital death” of 0.335, treating live discharge as censoring. Today’s estimate is 0.107. Using Lesko et al.’s language about measures of occurrence and competing events: which number describes the world these patients actually inhabit, and what hypothetical world does the other one describe?

YOUR ANSWER:

Q2. A collaborator proposes reporting the age-and-MAP-adjusted sex odds ratio (1.34) as “the” sex difference in mortality, because, noting that the AIC statistic for the null model (AIC = 867.47) is larger than the AIC statistic for the adjusted model (AIC = 793.62), “it’s the better model.” Is this reasoning correct? What question is the adjusted estimate answering? State one legitimate descriptive reason you might nevertheless prefer an adjusted (directly standardized or reweighted) estimate over the crude 1.21.

YOUR ANSWER:

Q3. Name two assumptions your descriptive claim about the target population needs, and one assumption it does not need. For the latter, say in one sentence why not.

YOUR ANSWER:

4 Step 2 — Reading it predictively (~11 min)

Underlying question:

“At the 60-minute landmark, can we flag the patients most likely to die in hospital, so the unit can watch them more closely?”

Vansteelandt and Steen call the target of a prediction analysis the predictimand, and it has its own required elements: the population in which predictions will be made (not the one you trained, or fit, your model in, but the next one), the moment of prediction (the predictors available at that moment you will make a prediction in the population, and the loss function (the measure that you will use to determine how “good” / “bad” a prediction is)

Whether we determine that a prediction is good or bad must come from the use of the prediction, and NOT from the outputs of a regression modeling software algorithm.

Your predictimand — complete the template

For ________________________________ (whom — in what deployment population?), at ________________________ (what moment?), predict ________________________ (what outcome) using ________________________________ (which inputs, available at that moment), judged by ________________________ (what measure of predictive performance, and why that one).

YOUR ANSWER:

One important check you must make before computing a performance number is whether every predictor you will use is knowable at the landmark when you will generate the prediction. Age, sex, baseline MAP, baseline heart rate, are all recorded at or before time zero. early_vasopressor is subtler: it is defined by behavior during (t0, t0+60], so it becomes fully known exactly at the landmark, and thus can be used.

Something like followup_days, however, would not be at the time of the landmark: it happens in the future. Thus, while it will likely greatly improve the predictive accuracy of the model, you simply cannot use it for the purposes of the predictimand you are targeting.

Now evaluate the model as a forecaster. Complete TODO (2a) — refit without the exposure. Everything else is provided, including a three-line AUC function: the AUC is the probability that a randomly chosen death was assigned a higher predicted risk than a randomly chosen survivor (0.5 = coin flip, 1 = perfect ranking).

# AUC = P(random death outranks random survivor in predicted risk)
auc_fun <- function(p, y) {
  r <- rank(p); n1 <- sum(y == 1); n0 <- sum(y == 0)
  (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

# Do consciously what glm() did silently in Step 0: form the analytical
# sample of 1,271 with a complete heart rate.
dat_model <- coh |> drop_na(baseline_hr)

p_full <- predict(fit, newdata = dat_model, type = "response")
c(mean_predicted = mean(p_full),
  observed_risk  = mean(dat_model$hospital_death))
auc_fun(p_full, dat_model$hospital_death)

# TODO (2a): refit WITHOUT early_vasopressor (same 1,271 rows, so the
# comparison is fair), and compute its AUC
fit_noexp <- glm(hospital_death ~ ???,
                 family = binomial, data = dat_model)
auc_fun(predict(fit_noexp, type = "response"), dat_model$hospital_death)

# and the exposure by itself, for reference
fit_only <- glm(hospital_death ~ early_vasopressor,
                family = binomial, data = dat_model)
auc_fun(predict(fit_only, type = "response"), dat_model$hospital_death)

Checkpoints: mean predicted risk 0.107 = observed risk 0.107; AUC of the full model 0.7246; without early_vasopressor 0.7225; early_vasopressor alone 0.5128.

First, note that early_vasopressor, the variable with the largest coefficient in the printout (0.429), is worth only 0.002 units of AUC. Alone, it ranks patients barely better than a coin flip. This is true even though the conditionally adjusted OR for early_vasopressor is quite large. “Importance” measured by coefficient size and importance for prediction are different quantities.

Second, that perfect calibration line (0.107 = 0.107) is not evidence of anything: a logistic model fit by maximum likelihood reproduces the overall rate in its own training data by construction. It is a check that cannot fail, and it is thus worthless. Instead, to make some meaning of this calibration check, you’d have to re-fit the model on a subset of the data (a training sample), and then predict in the hold-out (validation) sample, constructing your calibration measure in that hold-out sample.

Third, and importantly for prediction, nothing we computed is validated in an out of sample validation or hold-out set: every number above is in-sample. A real predictive analytic tool would require out-of-sample performance measures, in data resembling the deployment population. We have not done that here, and this should be clear.

Questions

Q4. Your colleague ranks the printout’s predictors by p-value and declares age and heart rate “the important predictors” and the rest “noise.” Give one reason why this is still the wrong procedure for a predictive question.

YOUR ANSWER:

Q5. Suppose the clinical unit deploys your predictive flags, and clinicians respond by starting vasopressors earlier and more aggressively in the flagged patients. Next year, the model is recalibrated against observed deaths. What would happen to the performance of the algorithm, and why is this a problem the training data cannot provide information on? (Vansteelandt and Steen: this is where a predictive question quietly acquires a causal ingredient.)

YOUR ANSWER:

Q6. The eICU demo is from 2014–2015. Name two concrete ways the joint distribution of (predictors, outcome) could differ in the population where this model would actually be deployed, and what would this do to the model’s usefulness.

YOUR ANSWER:

5 Step 3 — Reading it causally (~11 min)

In Lab 1, you wrote the causal contrast this dataset was built to support: the difference (or ratio) in in-hospital mortality risk had everyone followed early initiation versus had everyone followed no early initiation. These are two precisely worded strategies that can be implemented in the landmark population.

In potential-outcome notation, this contrast compares \(E(Y^{a=1})\) with \(E(Y^{a=0})\). But note that this estimand is model agnostic: there’s no mention of a model, a link function, or a coefficient, or any other related concept.

Can the printout’s early_vasopressor row (logOR 0.429) be interpreted as that contrast?

In the lecture, we briefly mentioned the need to separate causal questions into two separate “layers”, so to speak:

The identification layer belongs to the question and the design, not to the model. Here, we consider assumptions such as exchangeability (do age, sex, baseline MAP, and heart rate capture everything that made clinicians reach for norepinephrine within the hour?), positivity (27 exposed patients, 2.1% of the cohort), and consistency (is the strategy wording from Lab 1’s 6b precise enough that “the outcome under early initiation” means one thing?).

The model layer is added by our choice of estimator, on top of identification: that the exposure’s effect is a constant log-odds shift across every age–sex–MAP–HR combination, and that the confounders enter linearly on the log-odds scale. These constraints have nothing to do with the causal question; they come with the fitting routine we used.

Run the chunk (provided) and look at what the design layer costs.

# the only row of the printout the causal question is about
b  <- coef(fit)["early_vasopressor"]
se <- summary(fit)$coefficients["early_vasopressor", "Std. Error"]
round(c(logOR = b, SE = se, OR = exp(b),
        lcl = exp(b - 1.96 * se), ucl = exp(b + 1.96 * se)), 2)

# and the raw material behind it
with(coh, table(early_vasopressor, hospital_death))

Checkpoints: conditional OR 1.54, 95% CI 0.57 to 4.16, which is an interval compatible with “roughly halves the odds” to “quadruples them.” The 2×2 table is the one you built in Lab 2 Step 1: 6 deaths among 27 exposed, 130 among 1,247 unexposed (crude OR 2.45; Lab 2’s crude risk ratio was 2.13). The adjusted estimate is smaller than the crude one, and plagued with uncertainty.

Questions

Q7. Sort the following assumption statements into the identification layer, the model layer, or neither, if one belongs to a different conversation entirely:

(a) Clinicians’ reasons for initiating within the hour are fully captured by age, sex, baseline MAP, and heart rate. (b) The exposure shifts the log-odds of death by the same amount at every age, sex, MAP, and HR. (c) Patients of every covariate pattern in the landmark population had a real chance of both initiating early and not. (d) The log-odds of death is linear in age. (e) The Lab 1 strategy definitions are precise enough that “the outcome under early initiation” is well-defined. (f) The 18 patients with blank discharge status are missing at random.

YOUR ANSWER:

The takeaway. The coefficient table never really changed. What changed, three times, however, was the question (i.e., the estimand), the parts of the printout that mattered, and the assumptions we need to make to interpret things the way we want to.

Doing this in the wrong direction (“fit the model, then interpret the coefficients in a descriptive, predictive, or causal framework”) is commonplace.

Where this goes next: Week 5 dissects the model itself, including the left-hand side, the right-hand side, and the distinction between target and nuisance functions. And then week 6 builds up our understanding of generalized linear models, and what exactly happens behind the scenes when we use them.

References

Carlin JB, Moreno-Betancur M. On the Uses and Abuses of Regression Models: A Call for Reform of Statistical Practice and Teaching. Stat Med. 2025;44(13–14):e10244.
Vansteelandt S, Steen J. Discussion of “On the Uses and Abuses of Regression Models: A Call for Reform of Statistical Practice and Teaching.” Stat Med. 2025;44:e10312. (Included in the course PDF of the Carlin & Moreno-Betancur paper.)
Lesko CR, Fox MP, Edwards JK. A Framework for Descriptive Epidemiology. Am J Epidemiol. 2022;191(12):2063–2070.
Greenland S. Some Ways to Make Regression Modeling More Helpful Than Misleading. Stat Med. 2025;44:e10313.
Pollard T, et al. eICU Collaborative Research Database Demo (v2.0.1). PhysioNet, 2021. doi:10.13026/4mxk-na84 (the Lab 1 source data).