Lab 4 — Dissect a Regression Model

EPID 785R · Regression and Study Design

Author

YOUR NAME HERE

Published

October 4, 2026

Introduction

This the last lecture, we took apart a regression model: a left-hand side, a right-hand side that splits into a target and a nuisance function, a link, a distribution, and (sometimes) an offset, with the parametric-to-nonparametric spectrum describing how much structure the right-hand side imposes. Each of these is a choice, and each choice decides what the fitted numbers mean.

In this lab we explore some of the implications of turning these dials one at a time, on the same cohort, holding everything else fixed, and watch which numbers move and which do not. We spend some time figuring out what these changes mean. Like last week, the lab is short on new code and long on concept: most chunks are provided and most work is in the questions.

In this lab, you will:

  1. inventory the parts of the model you fit in Lab 3, by asking the fitted object for its own anatomy;
  2. change the link (identity → log → logit) while holding the left-hand side, the right-hand side, and the distribution fixed;
  3. add an offset to a Poisson model and watch a risk ratio become a rate ratio; verify the result against the hand calculation from Lab 2; and see what the offset literally is by letting its coefficient go free;
  4. fit one relationship (in-hospital death against age) three ways: linear, natural cubic spline, and saturated (the nonparametric maximum likelihood estimator), and check that the saturated fit is nothing more than the stratum proportions;
  5. visualize the spectrum these three fits span, read off what each position costs and buys, and watch what the choice does to a target coefficient sitting next to the age term.

Deliverables (yours to keep, nothing is submitted): this document rendered to HTML with your answers in the answer blocks. Four chunks contain TODOs with ??? placeholders (they error until fixed). When they all run, flip eval: false to eval: true in the YAML header and re-render.

Background for this lab: the Week 5 notes (The Anatomy of a Regression Model) and Clark & Berry’s Models Demystified, chapters 2–4 (reading list). Greenland’s (2025) commentary on Carlin & Moreno-Betancur supplies the lab’s organizing image, a fitted model as a lossy compression of the data.

1 Before you start (~3 min)

  1. Open lab4.Rproj in RStudio.
  2. The lab needs your Lab 1 analytic dataset. As in Labs 2 and 3, a copy ships in this lab’s data/ folder in the starter package.
  3. Packages: the tidyverse plus splines (which ships with base R; it provides the natural cubic spline in Step 3).
pacman::p_load(tidyverse, splines)

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

# As in Labs 2 and 3: drop the 18 stays with a blank discharge status
# (a complete-case analysis, for convenience, not as a recommendation).
coh <- lab1 |> filter(!is.na(hospital_death))
nrow(coh)    # expect 1,274

No model in this lab uses baseline_hr, so (unlike Lab 3) every fit below uses all 1,274 rows.

2 Step 0 — The parts list (~4 min)

Here is the model from Lab 3, Step 0, refit. Before changing any dials, name the regression model parts. A fitted glm object will tell you most of them (the chunk is provided):

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

formula(fit_lab3)            # the left- and right-hand sides
fit_lab3$family$family       # the distribution
fit_lab3$family$link         # the link
is.null(fit_lab3$offset)     # TRUE: no offset
length(coef(fit_lab3))       # how many parameters
NoteThe anatomy of Lab 3’s model — complete the table
part your answer
Left-hand side: the outcome, and the summary of it being modeled
Link: the scale on which the right-hand side is additive
Distribution, and the variance it implies for each observation
Right-hand side: the variables, and the functional form each enters with
Under Lab 3’s causal reading, the target part of the right-hand side
Under Lab 3’s causal reading, the nuisance part of the right-hand side
Offset
What exp(coef) for early_vasopressor is, as a contrast (a difference? a ratio? of what, comparing whom?)

The rest of the lab will hold the left-hand side fixed (in-hospital death, as a probability), hold the distribution fixed within each step, and change one dial at a time: the link (Step 1), the offset (Step 2), and the flexibility of the right-hand side (Steps 3 and 4).

4 Step 2 — Add an offset (~15 min)

The lecture’s definition: an offset is a term on the right-hand side whose coefficient is fixed at one rather than estimated. Its canonical use is to turn a log-link count model into a rate model, because on the log scale, subtracting \(\log T\) from both sides divides the mean by \(T\):

\[\log E(Y \mid X, T) = \beta_0 + \beta_1 X + \log T \quad\Longleftrightarrow\quad \log \frac{E(Y \mid X, T)}{T} = \beta_0 + \beta_1 X.\]

The cohort has everything a rate model needs: an event count per person (hospital_death, which is 0 or 1) and that person’s time at risk (followup_days, from the landmark to hospital discharge, dead or alive).

4.1 Without the offset: Lab 2’s risk ratio

Begin with a Poisson model with a log link and no offset (provided). You fit exactly this in Lab 2 to get the full-cohort risk ratio, where its model-based standard error was the one we replaced with a robust one.

fit_pois <- glm(hospital_death ~ early_vasopressor,
                family = poisson(link = "log"), data = coh)
round(exp(coef(fit_pois)), 4)   # the unexposed risk, and the risk ratio

Checkpoint: exp(intercept) = 0.1043, the risk among the unexposed (130/1,247), and exp(coefficient) = 2.13, Lab 2’s risk ratio. Compare with Step 1: the Poisson distribution with a log link gives the same risk ratio that the binomial distribution with a log link gave (2.13), and, since the model is saturated, the same fitted risks. The distribution dial and the link dial really are separate.

4.2 By hand: rates and the rate ratio

Now compute what a rate model targets, by hand, as you did in Lab 2. Complete TODO (2a).

# TODO (2a): deaths and person-days by exposure group, then rates
rates <- coh |>
  group_by(early_vasopressor) |>
  summarise(n      = n(),
            deaths = sum(???),
            pdays  = sum(???),
            risk   = deaths / n,
            rate   = deaths / pdays)        # deaths per person-day
as.data.frame(rates)
round(100 * rates$rate, 2)                  # per 100 person-days, as in Lab 2
c(RR = rates$risk[2] / rates$risk[1], IRR = rates$rate[2] / rates$rate[1])

Checkpoints: 130 deaths over 7,310.1 person-days among the unexposed and 6 over 168.3 among the exposed (7,478.4 in total, Lab 2’s number): rates of 1.78 and 3.56 per 100 person-days, a rate ratio of 2.00. The risk ratio is 2.13.

4.3 With the offset: the same machine, now a rate model

Add the offset. Complete TODO (2b): the term is offset(log(followup_days)), and nothing else in the call changes.

# TODO (2b): the same Poisson model, with log person-time as an offset
fit_rate <- glm(hospital_death ~ early_vasopressor + offset(???),
                family = poisson(link = "log"), data = coh)
round(coef(fit_rate), 4)
round(exp(coef(fit_rate)), 4)           # the unexposed rate (per day), and the rate ratio
round(100 * exp(coef(fit_rate)[1]), 3)   # the unexposed rate per 100 person-days

Checkpoints: exp(intercept) = 0.01778 deaths per person-day (1.78 per 100 person-days, the unexposed rate you just computed) and exp(coefficient) = 2.004: the hand-computed rate ratio, exactly. (The coefficient’s model-based standard error, 0.4175, is practically unchanged from the no-offset model’s 0.4176.)

4.4 What the offset literally is

An offset is a covariate whose coefficient is fixed at one. To see that, remove the fixing: put log(followup_days) on the right-hand side as an ordinary covariate and let the data estimate its coefficient (provided).

fit_free <- glm(hospital_death ~ early_vasopressor + log(followup_days),
                family = poisson(link = "log"), data = coh)
round(coef(fit_free), 4)
round(exp(coef(fit_free)["early_vasopressor"]), 3)
c(deviance_offset = deviance(fit_rate), deviance_free = deviance(fit_free))

# for the record: the follow-up each group accumulated, by outcome
coh |>
  group_by(early_vasopressor, hospital_death) |>
  summarise(n = n(), median_days = round(median(followup_days), 2), .groups = "drop")

Checkpoints: set free, the coefficient on log follow-up is −0.235, not 1, and the exposure coefficient moves to 0.639 (exponentiated: 1.89); the deviance drops, from 863.0 to 592.9. The table says why the sign is negative: in this cohort follow-up ends at death or discharge, so the patients who die are the ones with short follow-up (median 2.75 days among unexposed decedents vs. 3.92 among unexposed survivors; 1.73 vs. 7.94 among the exposed). When we include it in the model as a covariate (meaning, we estimate a parameter for it) the time variable stops measuring time at risk and starts acting as a proxy for the outcome. Because it is so highly predictive of the outcome, it absorbs much of the outcome’s variation in the model.

The rate model asserts instead that expected deaths scale in proportion to time at risk.

One last check of what the offset model actually uses (provided): collapse the cohort to two rows, total deaths and total person-days per group, and refit.

two_rows <- coh |>
  group_by(early_vasopressor) |>
  summarise(deaths = sum(hospital_death), pdays = sum(followup_days))
as.data.frame(two_rows)

fit_rate_2 <- glm(deaths ~ early_vasopressor + offset(log(pdays)),
                  family = poisson(link = "log"), data = two_rows)
round(rbind(individual = coef(fit_rate), aggregated = coef(fit_rate_2)), 6)

Checkpoint: the two-row fit reproduces the 1,274-row coefficients to every printed decimal. For this model the individual rows carry no information beyond the group totals of deaths and person-time: the rate ratio is the ratio of those totals.

NoteQuestions

Q4. The no-offset model returned the risk ratio (2.13) and the offset model the rate ratio (2.00). For each, write one sentence naming the denominator it uses. Then use the median follow-up table to explain why the rate ratio is the smaller of the two in this cohort.

YOUR ANSWER:

Q5. When unconstrained (i.e., not fixed at 1), the coefficient on log follow-up time was −0.23, and the deviance fell from 863 to 593. Why would letting the data choose this coefficient be a mistake in a rate model, even though the fit is so much “better”? What assumption does fixing it at 1 add, and is that assumption plausible over a hospital stay?

YOUR ANSWER:

Q6. The two-row dataset reproduced the full fit exactly. Suppose age (in years) were added to the rate model. Would two rows still suffice? Describe the smallest aggregated dataset that would reproduce that fit, and say in one sentence why.

YOUR ANSWER:

5 Step 3 — One relationship in three ways (~15 min)

The right-hand side’s last dial is its flexibility: how lossy the compression is, that is, how many parameters the model spends and therefore how much of the pattern between the outcome and a covariate it keeps. The lecture located three positions on a spectrum, parametric, semiparametric, and nonparametric, and showed that a saturated model (one parameter per covariate pattern) is the nonparametric end: a lossless compression that keeps every covariate pattern’s proportion and filters out nothing, whose maximum likelihood estimator is called the NPMLE. Step 1 varied which pattern a two-parameter summary keeps; this step varies how much is kept.

Age lets us see the whole spectrum in one relationship: 73 distinct values, with between 2 and 59 patients at each age value. This is enough to fit a model with a parameter for every unique year of age.

So we’ll fit in-hospital death against age three ways, all with the logit link (so “linear” means linear in age on the log-odds scale; the figure in Step 4 translates everything back to the risk scale). The chunk is provided; its three right-hand sides are:

  • linear: age enters as itself (one slope);
  • natural cubic spline with 4 degrees of freedom, ns(age, df = 4) from the splines package: a smooth curve made of cubic pieces joined at three interior knots (placed at the age quartiles, 55, 68, and 78) and constrained to be linear beyond the data’s ends;
  • saturated: age as a factor, factor(age), one parameter per observed year of age.
# the same relationship, three ways (provided)
fit_lin <- glm(hospital_death ~ age,              family = binomial, data = coh)
fit_spl <- glm(hospital_death ~ ns(age, df = 4),  family = binomial, data = coh)
fit_sat <- glm(hospital_death ~ factor(age),      family = binomial, data = coh)

spectrum <- data.frame(
  model      = c("linear", "natural spline (df = 4)", "saturated"),
  parameters = c(length(coef(fit_lin)), length(coef(fit_spl)), length(coef(fit_sat))),
  deviance   = round(c(deviance(fit_lin), deviance(fit_spl), deviance(fit_sat)), 2),
  AIC        = round(c(AIC(fit_lin), AIC(fit_spl), AIC(fit_sat)), 2))
spectrum

Checkpoints: 2, 5, and 73 parameters; residual deviances 820.60, 815.95, 760.23; AICs 824.60, 825.95, 906.23. (The null deviance, for reference, is 865.47.) Deviance falls at every step, as it must: each model contains the previous one. AIC, which charges two units per parameter, rises at every step.

Now let’s check the lecture’s claim that the saturated fit is the NPMLE and reproduces the data. Check it (provided): the saturated model’s fitted risk at each age should be exactly the observed proportion of deaths at that age.

by_age <- coh |>
  group_by(age) |>
  summarise(n = n(), deaths = sum(hospital_death), p_obs = mean(hospital_death))

by_age$p_sat <- predict(fit_sat, newdata = data.frame(age = by_age$age),
                        type = "response")
max(abs(by_age$p_obs - by_age$p_sat))        # zero, to numerical precision

# the same saturated model under the identity link, fit by least squares
by_age$p_sat_lm <- predict(lm(hospital_death ~ factor(age), data = coh),
                           newdata = data.frame(age = by_age$age))
max(abs(by_age$p_obs - by_age$p_sat_lm))     # also zero

# where the saturated model is "estimating" from very little
by_age |> filter(age %in% c(18, 28, 44, 85, 88, 90))

# the saturated printout, for the ages with no deaths
round(summary(fit_sat)$coefficients[c("(Intercept)", "factor(age)19", "factor(age)28"), ], 2)

Checkpoints: both differences are zero to numerical precision (about 1e−8 for the logistic fit, which stops iterating when it is close enough, and about 1e−14 for least squares). The saturated model is the stratum-specific proportions, and (echoing Step 1) it does not matter which link you fit it with. The filtered rows show how those proportions are calculated: 0 of 7 at age 18, 1 of 10 at 28, 0 of 3 at 44, 6 of 15 at 85, 1 of 13 at 88, and 14 of 59 at 90 (which is where Lab 1 placed everyone over 89).

Also, look at the last printout. Thirty of the 73 ages contain no deaths; their fitted risk is 0, whose log-odds is −∞, and the printout shows it as an intercept near −18.6 (age 18 is the reference level) with a standard error in the thousands. Age 28, with one death, gets a coefficient of +16.4 to climb back. In next week’s lecture under the name separation, we’ll see why this is produced here (not by an unusual covariate but by asking for one parameter per age from a cohort that cannot supply it).

NoteQuestions

Q7. Deviance falls from the linear to the spline to the saturated model, and AIC rises. In your own words, what is each model “paying” with, and what is it “buying”? Then: the spline’s three extra parameters reduced the deviance by 4.65 (a likelihood-ratio p-value of about 0.20). Does that settle whether the spline is the better model? Better for what?

YOUR ANSWER:

Q8. The saturated fit says mortality is 0.40 at age 85 and 0.08 at age 88. Is the five-fold drop a feature of the cohort or of the estimator? If you wanted to move away from the saturated end just far enough to repair it, name a model that sits between the saturated fit and the spline, and say what it filters out that the saturated fit kept.

YOUR ANSWER:

6 Step 4 — Visualize the spectrum (~12 min)

Draw the three fits over the observed proportions (provided; point size is the number of patients at each age).

grid    <- data.frame(age = seq(18, 90, by = 0.5))
sat_grd <- data.frame(age = sort(unique(coh$age)))   # a factor can only be predicted at its levels

spectrum_fits <- bind_rows(
  tibble(age = grid$age,    risk = predict(fit_lin, grid,    type = "response"),
         Model = "Linear (2 parameters)"),
  tibble(age = grid$age,    risk = predict(fit_spl, grid,    type = "response"),
         Model = "Natural spline (5 parameters)"),
  tibble(age = sat_grd$age, risk = predict(fit_sat, sat_grd, type = "response"),
         Model = "Saturated (73 parameters)")) |>
  mutate(Model = factor(Model, c("Linear (2 parameters)",
                                 "Natural spline (5 parameters)",
                                 "Saturated (73 parameters)")))

ggplot() +
  geom_point(data = by_age, aes(x = age, y = p_obs, size = n),
             color = "grey60", alpha = .7) +
  geom_line(data = spectrum_fits, aes(x = age, y = risk, color = Model),
            linewidth = .8) +
  scale_color_manual(values = c("#009E73", "#D55E00", "#0072B2")) +
  scale_size_area(max_size = 5, guide = "none") +
  labs(x = "Age (years)", y = "Risk of in-hospital death", color = NULL) +
  theme_classic() + theme(legend.position = "bottom")

Evaluate this figure against the lecture’s kidney figure. The linear fit (linear on the logit scale, so gently curved here) is the most compressed of the three and filters out whatever its form cannot carry: it cannot see the flat stretch below age 40 or the steepening after 70, both of which live in its residuals. The spline follows both, and borrows strength across neighboring ages to do it. The saturated fit passes through every observed proportion, including the thirty zeros and the 0.40 at age 85 that rests on 15 patients. It is a lossless compression of the data, in that it compresses nothing.

6.1 The target next to the nuisance

Finally, let’s bring the exposure back. Under Lab 3’s causal reading, the exposure coefficient is the target function and age is a nuisance function: we carry it only so the exposure contrast is adjusted for it. Here is what the three positions on the spectrum do to the target when age sits beside it (provided). The adjustment set is deliberately age alone, so that the only thing changing is age’s functional form; Lab 3’s fuller model, with sex, MAP, and heart rate as well, gave 1.54.

or_row <- function(fit) {
  b  <- unname(coef(fit)["early_vasopressor"])
  se <- unname(sqrt(vcov(fit)["early_vasopressor", "early_vasopressor"]))
  c(parameters = length(coef(fit)), logOR = b, SE = se, OR = exp(b),
    lcl = exp(b - 1.96 * se), ucl = exp(b + 1.96 * se))
}
round(rbind(
  age_linear    = or_row(glm(hospital_death ~ early_vasopressor + age,
                             family = binomial, data = coh)),
  age_spline    = or_row(glm(hospital_death ~ early_vasopressor + ns(age, df = 4),
                             family = binomial, data = coh)),
  age_saturated = or_row(glm(hospital_death ~ early_vasopressor + factor(age),
                             family = binomial, data = coh))), 3)

# how many of the 73 age strata contain an exposed patient at all?
coh |>
  group_by(age) |>
  summarise(exposed = sum(early_vasopressor)) |>
  summarise(strata_with_an_exposed_patient = sum(exposed > 0))

Checkpoints: OR 2.16 (SE 0.483) with age linear, 2.14 (SE 0.482) with the spline, 2.37 (SE 0.517) with age saturated; only 19 of the 73 age strata contain an exposed patient.

NoteQuestions

Q9. Which of the three curves is “right”? Answer three times, once for each of Lab 3’s questions: (i) a descriptive question whose descriptimand is “the proportion dying, by single year of age, in this cohort”; (ii) a predictive question about a patient who will arrive next year at age 61; (iii) the causal question, where age is a nuisance variable. For each, say which position on the spectrum you would choose and what you would be assuming.

YOUR ANSWER:

Q10. Going from a linear to a spline age term moved the exposure OR from 2.16 to 2.14 and left its SE at 0.48; going to a saturated age term moved it to 2.37 and the SE to 0.52. Which of these moves is the nuisance function buying bias protection, and which is it paying variance? Why does saturation hurt here in particular? (Hint: the last line of the chunk.)

YOUR ANSWER:

The takeaway. Nothing about the cohort changed today. We turned three dials of one model and watched. The link, which in a lossless (saturated) model only chooses the coordinates in which the same two risks are stored, began deciding which feature of the data survives the moment the compression became lossy. The offset, a coefficient we fix rather than estimate, moved the same Poisson machine from a risk ratio to a rate ratio, and would have been “estimated” to mean something else entirely had we let it go free. And the flexibility of the right-hand side, where every step toward the nonparametric end made the compression less lossy, lowered the deviance, and raised the price in parameters, until the model stopped summarizing and started repeating the data. Each setting is a choice, none of them is in the data, and each one decides what the number you report means.

Where this goes next: Tuesday’s lecture turns the fourth dial, the distribution, and opens the fitting machine (least squares, maximum likelihood, and the iteratively reweighted algorithm behind every glm() call you ran today). Lab 5 fits the whole GLM menu to this same outcome, including the identity- and log-link failures previewed in Step 1, and adds the diagnostics.

References

  • Greenland S. Some Ways to Make Regression Modeling More Helpful Than Misleading. Stat Med. 2025;44:e10313. (The lossy-compression image; in the course PDF of Carlin & Moreno-Betancur.)
  • Clark M, Berry S. Models Demystified: A Practical Guide from Linear Regression to Deep Learning, chapters 2–4. https://m-clark.github.io/book-of-models/
  • Naimi AI, Whitcomb BW. Estimating Risk Ratios and Risk Differences Using Regression. Am J Epidemiol. 2020;189(6):508–510.
  • Zou G. A Modified Poisson Regression Approach to Prospective Studies with Binary Data. Am J Epidemiol. 2004;159(7):702–706.
  • Greenland S, Mansournia MA, Altman DG. Sparse Data Bias: A Problem Hiding in Plain Sight. BMJ. 2016;352:i1981. (Background for Q10.)
  • Pollard T, et al. eICU Collaborative Research Database Demo (v2.0.1). PhysioNet, 2021. doi:10.13026/4mxk-na84 (the Lab 1 source data).