Lab 4 — Instructor answer key

Companion file: solution.R (in this folder), the lab chunk code with the four TODOs (1a, 1b, 2a, 2b) completed and the Step 3 chunk as the handout now provides it, in chunk order. Run it from the lab root (Rscript solutions/solution.R) and it prints every checkpoint in the handout plus the extras quoted below (confidence intervals, likelihood ratio tests, the decade-bin model, fitted risks at age 61, and the identity-link boundary fit).

Contents: TODO code answers · the parts list (Step 0) · expected checkpoints · model answers Q1–Q10 · understanding notes

TODO code answers

TODO	Answer
1a	As edited by the instructor on 2026-10-04 the link names are given; the blanks are `coef(fit_id)`, `coef(fit_log)`, `coef(fit_lgt)` in the `rbind()`, then `exp(coef(fit_log))` and `exp(coef(fit_lgt))`, then `fitted(fit_id)`, `fitted(fit_log)`, `fitted(fit_lgt)` inside the `sort(unique(round(...)))` calls. (The `# TODO (1a)` comment and the “filling in the other two link names” sentence still describe the old blanks.) Nothing else changes: same outcome, same right-hand side, same distribution. The three fits are saturated (two parameters, two groups), so their fitted values and deviances are identical and only the coefficient’s meaning differs: RD 0.118, log RR 0.757 (RR 2.13), log OR 0.898 (OR 2.45).
1b	`hospital_death ~ age` under both links. The two fits converge in 6 and 5 IRLS iterations. The provided identity-link call is meant to fail: with zero deaths at ages 18–27, the straight line’s MLE lies on the boundary of the parameter space (risk exactly 0 at age 18), and `glm()` cannot find valid starting values. Students who “fix” it by supplying `start =` and a large `maxit` (e.g. `start = c(0.001, 0.0017), control = glm.control(maxit = 500)`) get a converged boundary solution, slope 0.00229 and fitted risk 0.000000 at age 18. That is a fine aside, not the lesson.
2a	`deaths = sum(hospital_death)`, `pdays = sum(followup_days)`. Lab 2’s Step 1 computed the same quantities with `a1 / pt1` notation; the rates are 1.78 and 3.56 per 100 person-days and the rate ratio 2.004.
2b	`offset(log(followup_days))`. Equivalent spellings: the `offset = log(followup_days)` argument to `glm()`, or `+ log(followup_days)` with its coefficient constrained — which is what the offset does. The exponentiated coefficient equals the hand rate ratio to every printed decimal.
3a	Provided in the handout (the instructor filled the three blanks on 2026-10-04, and the TODO labels were removed the same day). If ever re-blanked, the answers are `hospital_death ~ age`, `~ ns(age, df = 4)`, `~ factor(age)`. The spline’s three interior knots sit at the age quartiles (55, 68, 78) with boundary knots at 18 and 90; `ns(age, df = 4)` contributes 4 columns, hence 5 parameters with the intercept. `factor(age)` contributes 72 columns (73 levels, age 18 as reference).

The parts list (Step 0)

Model answers for the table; any phrasing that names the right object is fine.

part	model answer
Left-hand side	`hospital_death`, coded 1/0; the summary modeled is its conditional mean, which for a binary outcome is the conditional probability (risk) of in-hospital death given the covariates.
Link	logit: the right-hand side is additive on the log-odds scale, $\log[P/(1-P)]$.
Distribution	binomial (Bernoulli, one trial per row); implied variance $\mu_i(1-\mu_i)$ for each patient, no free variance parameter.
Right-hand side	an intercept; `early_vasopressor` (0/1 indicator); `age` (linear, per year); `sex` (indicator, Female reference); `baseline_map` and `baseline_hr` (linear, per unit); no product terms, so every term is additive on the log-odds scale. 6 parameters.
Target (causal reading)	the coefficient on `early_vasopressor`, $\alpha_1$.
Nuisance (causal reading)	everything else: the intercept and the age, sex, MAP, and HR terms, the function $\eta(\mathbf{C})$ carried so that the exposure contrast is adjusted.
Offset	none (`fit_lab3$offset` is `NULL`).
`exp(coef)` for the exposure	a conditional odds ratio: the odds of in-hospital death among early initiators divided by the odds among non-initiators of the same age, sex, baseline MAP, and baseline HR, assumed constant across those strata (1.54 in Lab 3). A ratio, not a difference; of odds, not risks; conditional, not marginal.

Expected checkpoints

Verified 2026-10-04 against the canonical Lab 1 dataset (1,292 rows → 1,274 complete-case; no model in this lab drops further rows). No seeds, no randomness anywhere.

setup:   nrow(lab1) = 1292; nrow(coh) = 1274
Step 0:  formula hospital_death ~ early_vasopressor + age + sex + baseline_map + baseline_hr
         family binomial, link logit, offset NULL, 6 parameters
Step 1a: exposure-only, binomial, three links
         link      intercept   early_vaso   exp(early_vaso)   SE(early_vaso)
         identity    0.1043      0.1180        —  (RD 0.118)     0.080
         log        -2.2610      0.7569       2.132 (RR)         0.369
         logit      -2.1509      0.8981       2.455 (OR)         0.472
         fitted values 0.1042 / 0.2222 under all three (= 130/1247, 6/27)
         deviance 862.40 under all three; inverse links by hand all 0.2222
         aside box (delta method, logistic vcov -> log RR -> RD): log RR SE 0.3695 = glm 0.3695;
           RD SE 0.0805 = glm 0.0805. Key-only extras: hand formulas 0.0805 / 0.3695 /
           0.4721 (Woolf); glm logOR SE 0.472093 vs Woolf 0.472095 (last-IRLS-iteration weights)
Step 1b: death ~ age
         log:   intercept -4.7726, age 0.0363 -> RR 1.037/yr, 1.44/decade (6 iterations)
         logit: intercept -4.9599, age 0.0409 -> OR 1.042/yr, 1.51/decade (5 iterations)
         identity: "no valid set of coefficients has been found: please supply starting values"
         least squares: intercept -0.0894, slope 0.0030; 73 patients with fitted risk < 0
         fitted at 18/40/65/80/90: log 0.016/0.036/0.090/0.155/0.223
                                   logit 0.014/0.035/0.091/0.156/0.218
                                   LPM  -0.035/0.032/0.107/0.153/0.183
         AIC log 824.84, logit 824.60; 0 deaths among the 52 patients aged 18-27
Step 2:  Poisson no offset: exp(coef) 0.1043, 2.1316 (Lab 2's RR); SE 0.4176
         by hand: unexposed 130 / 7,310.06 pd = 1.78 per 100 pd
                  exposed     6 /   168.35 pd = 3.56 per 100 pd; total 7,478.4
                  RR 2.1316, IRR 2.0041
         offset model: coef -4.0295, 0.6952; exp 0.01778, 2.0041; SE 0.4175
                       95% CI for the IRR 0.88-4.54
         log(followup_days) free: coef -2.0139, 0.6387, -0.2348; exp(exposure) 1.894
                       deviance 863.03 (offset) vs 592.88 (free)
         median follow-up (days): unexposed alive 3.92 / dead 2.75; exposed alive 7.94 / dead 1.73
         mean follow-up: unexposed 5.862, exposed 6.235; IRR = RR x 5.862/6.235 = 2.0041
         two-row refit: coefficients identical to 1e-8 (SEs differ in the 5th decimal, see notes)
Step 3:  linear 2 params, deviance 820.60, AIC 824.60
         ns(age, 4) 5 params, deviance 815.95, AIC 825.95   (knots 55/68/78)
         factor(age) 73 params, deviance 760.23, AIC 906.23  (17 iterations, converged)
         null deviance 865.47
         LRT spline vs linear 4.65 on 3 df, p = 0.20; saturated vs spline 55.72 on 68 df, p = 0.86
         max |p_obs - saturated fit| 8.6e-9 (logit) / 2.6e-14 (least squares)
         by age: 18 -> 0/7; 28 -> 1/10; 44 -> 0/3; 85 -> 6/15; 88 -> 1/13; 90 -> 14/59
         30 ages (254 patients) have zero deaths; n per age min 2, median 15, max 59
         saturated printout: intercept -18.57 (SE 2465); factor(age)28 +16.37 (SE 2465)
         extras: decade bins (<30, 30s, ..., 80+) 7 params, deviance 814.51, AIC 828.51
                 fitted risk at age 61: linear 0.078, spline 0.080, saturated 0.077 (1/13); saturated at 60 = 0 (0/22)
                 exact 95% CI: age 85 (6/15) 0.16-0.68; age 88 (1/13) 0.00-0.36
Step 4:  exposure OR | age linear 2.163 (SE 0.483, CI 0.84-5.57, 3 params)
                     | age spline 2.137 (SE 0.482, CI 0.83-5.50, 6 params)
                     | age saturated 2.375 (SE 0.517, CI 0.86-6.54, 74 params)
         19 of 73 age strata contain an exposed patient (27 exposed, ages 36-90)

Model answers, Q1–Q10

Q1. The two models, in compression terms. The exposure-only model is lossless: two observed proportions, two parameters, so every link stores both numbers exactly and discards nothing. The link only chooses the coordinates (difference, ratio, odds ratio) in which the same pair is written, and decompressing with the inverse link returns the data unchanged, which is why fitted values and deviance are identical. With age, 73 proportions must be stored in two numbers, so the compression is lossy and something is filtered into the residuals. Now the link decides what is kept: the identity link keeps the straight-line component of risk against age (one risk difference per year), the log link the exponential component (one risk ratio per year), the logit the sigmoid component (one odds ratio per year). Three filters tuned to three patterns pass different things through, so the fitted values differ. The dividing line is saturation: once there are more covariate patterns than parameters, the link stops being a choice of coordinates and becomes a choice of what to discard. A student who says “the model stopped being saturated, so the link became a restriction on the shape of the curve” has the mechanism; the compression vocabulary (the lecture’s, Greenland’s, and this lab’s) adds the consequence, which is that the fitted model is a summary whose adequacy depends on which features the application needs, not on fit.

Q2. “Canonical” does not mean “required.” The canonical link is the one that makes the linear predictor equal the exponential family’s natural parameter (next week: it gives sufficient statistics, a concave log-likelihood, and IRLS = Newton’s method). It is a mathematical convenience of the distribution, not a requirement of the outcome. Step 1a shows three binomial fits with three links, all legitimate, all fitting equally well, each returning a different and correct contrast. The link is chosen for the scale on which you want the answer (RD, RR, OR), which comes from the question, not from the outcome being 0/1. In Greenland’s terms, a model form is adequate “to the extent it allows through those patterns” the application needs; which pattern to keep (a difference, a ratio, an odds ratio) is decided by the question, and the canonical link is merely the one whose arithmetic is most convenient for the fitting algorithm. Practical reasons to default to the logit anyway: its fitted values are confined to (0, 1), so the likelihood is always defined and the fit converges (Step 1b’s identity failure is exactly the alternative); and the odds ratio is estimable under outcome- dependent sampling (Week 3 / Week 11) where risks are not.

Q3. Two ways to a risk difference for age. (i) Use a model that did converge (log or logit) and take differences of fitted risks: e.g., risk at 66 minus risk at 65, or at any two ages of interest. This assumes the log- or logit-linear shape, under which the risk difference per year is not constant: about 0.0006 per year at age 20 and about 0.008 per year at age 85 under the logit fit. So “the RD per year” is several numbers, reported at chosen ages, or one number if averaged over the cohort’s age distribution (a marginal contrast, Week 8). (ii) Fit the straight line by least squares (the linear probability model in the figure) and report its slope, 0.0030 per year, with a robust standard error (as Lab 2 did for the exposure RD). This assumes risk is linear in age over the whole range, which the figure shows is false below age 30 (73 negative fitted risks). Also acceptable: the binomial identity-link fit with starting values (a boundary solution with risk 0 at age 18, slope 0.00229), restricting the age range to where a line is plausible, or an identity link with a spline in age (then the RD is again age-specific). Any answer should notice that the choice of link (Step 1) and the choice of functional form (Step 3) are both in play.

Q4. Risk ratio vs. rate ratio: denominators, and why 2.00 < 2.13. The risk ratio divides deaths by people: (6/27) / (130/1,247). The rate ratio divides deaths by person-days: (6/168.3) / (130/7,310.1). The two differ only through follow-up per person, and the identity is exact: IRR = RR × (mean follow-up, unexposed) / (mean follow-up, exposed) = 2.132 × 5.862/6.235 = 2.004. The exposed accumulated 6% more follow-up per person than the unexposed, which dilutes their rate relative to their risk. The median table explains where that time came from: exposed decedents died fast (median 1.73 days), but exposed survivors stayed long (median 7.94 days vs. 3.92 among unexposed survivors), and survivors are 21 of the 27 exposed, so they dominate exposed person-time. The two ratios answer different questions (Week 2’s risk vs. rate), and the offset is the dial between them.

Q5. Why not let the data choose the offset’s coefficient? Because with the coefficient unconstrained (estimated rather than fixed at 1) the model is no longer a rate model: it describes how the count of deaths varies with follow-up time, and in this cohort follow-up time is partly determined by the outcome (death ends it). The negative coefficient says “patients with more follow-up died less,” which is true and backwards: they had more follow-up because they did not die. The deviance falls from 863 to 593 because log follow-up is an excellent predictor of death, for that reason. Fitting the joint pattern of (death, time) better is not the question. Fixing the coefficient at 1 asserts that expected deaths are proportional to time at risk within a covariate pattern, i.e., a constant hazard over the stay (exponential survival within group), a relationship the model asserts rather than estimates. Plausibility: over a hospital stay after a hypotensive episode the hazard is probably highest in the first days and falls, so the assumption is at best approximate; the IRR is then a person-time-weighted average of time-varying rate ratios. Weeks 12–13 replace the constant-hazard assumption with survival models.

Q6. Would two rows still suffice with age in the model? No. The Poisson log-likelihood with an offset depends on the data only through the total events and total person-time within each distinct covariate pattern of the model (rows with the same linear predictor can be pooled without changing the likelihood). With exposure and age, the smallest aggregated dataset has one row per observed (exposure, age) combination, with its deaths and person-days: 73 unexposed rows plus the 19 ages that contain an exposed patient, 92 rows. Two rows no longer suffice because the model’s linear predictor now differs across ages, and the aggregation must respect every pattern the model distinguishes. (A sharper answer: with age entering linearly, the sufficient statistics are even smaller, the sums of deaths, deaths × exposure, and deaths × age, plus the person-time structure; but the glm call needs the cell-level rows to reproduce the fit.)

Q7. What each model pays and buys; does p = 0.20 settle it? Each step pays in parameters (degrees of freedom, and with them variance and the risk of fitting noise) and buys deviance, that is, fidelity to this sample. AIC prices the trade at two deviance units per parameter as an approximation to out-of-sample loss: the spline’s three extra parameters bought 4.65 of deviance (not quite the 6 they cost), and the saturated fit’s 68 extra parameters bought 55.7 (far short of 136). The LRT p-value of 0.20 says the curvature is not strongly demanded for in-sample fit. It settles nothing about which model is better, because “better” depends on the use: for the causal question the spline’s job is to protect the exposure coefficient from residual confounding by age, and Step 4 shows it does so at no cost in the target’s precision (SE 0.483 → 0.482), so it is the better model there even though AIC prefers the line; for a predictive question the criterion is out-of-sample performance, which neither the deviance nor the LRT measures; for a descriptive question by single year of age, neither model is the estimand. AIC and LRTs rank models on their fit to these data, not on their fitness for the question.

Q8. 0.40 at 85 vs. 0.08 at 88: cohort or estimator? The estimator. Six of 15 and one of 13 are proportions estimated from 15 and 13 patients; their exact 95% intervals are 0.16–0.68 and 0.00–0.36 and overlap substantially. Nothing in human physiology changes five-fold between 85 and 88, and the saturated model has no way to say so: it imposes no structure between neighboring ages, so it cannot borrow the 85-year-olds’ neighbors to steady the estimate. A model between the saturated fit and the spline: age in bins (a step function). Decade bins give 7 parameters, deviance 814.5, AIC 828.5. What it filters out is the year-to-year variation within a bin (the ages in a bin are pooled into one stored proportion); what it keeps, like the saturated fit, is any pattern across bins, since nothing constrains neighboring bins to resemble each other. Also acceptable: a saturated fit with shrinkage of the age effects toward a smooth curve (Week 10’s penalization, or a random-effects model), or a spline with more knots. The move buys stability at the price of the within-bin detail it discards, which is the bias–variance tradeoff of Week 4 with named positions.

Q9. Which curve is “right”? Three answers. (i) Descriptive, descriptimand = proportion dying by single year of age in this cohort: the saturated model is not an estimate of the descriptimand, it is the descriptimand, the by-age table under another name; no modeling assumption is needed (whether it describes a target population needs the sampling and measurement assumptions of Lab 3, Q3, and whether a single-year summary is useful depends on reporting it with its uncertainty, or coarsening it, which is a precision question, not a bias one). (ii) Predictive, a patient arriving next year at 61: the saturated fit predicts 0.077 (1 of 13), and would have predicted 0 for a 60-year-old (0 of 22), which is noise masquerading as precision; the linear and spline fits predict 0.078 and 0.080 and borrow strength from neighboring ages. Choose by out-of-sample loss (cross-validation, or a held-out year), and the assumption is that the smooth age–risk relationship transports to the deployment population (Lab 3’s predictimand discipline). (iii) Causal, age as nuisance: the aim is to remove confounding by age without spending precision on it, so a flexible-but-cheap nuisance form is the default: the spline. Assumed: that the exposure effect is constant across age on the logit scale (no product term), and that the spline absorbs enough of the age–outcome relationship that no residual confounding leaks into the exposure coefficient. The saturated nuisance is “assumption-free” about age’s form but, as Step 4 shows, costs precision and introduces sparse-data bias.

Q10. Bias protection vs. variance cost; why saturation hurts here. Linear → spline (2.16 → 2.14, SE unchanged at 0.48): this is the nuisance function buying bias protection, which in these data turned out to be nearly unnecessary (the age relationship is close to log-linear), and it cost nothing in precision. Spline → saturated (2.14 → 2.37, SE 0.48 → 0.52): the SE increase is the variance cost. The estimate moving away from the null is not bias protection but sparse-data bias (Greenland, Mansournia & Altman 2016): with one free intercept per age, the exposure contrast is estimated only inside the 19 strata that contain an exposed patient (strata with no exposed contribute nothing to the exposure score), most of those strata hold a single exposed patient, and unconditional logistic regression with many sparse strata is biased away from the null. Saturation hurts here in particular because the exposure is rare (27 patients) and scattered across 19 of 73 strata: the other 54 strata’s information about age is irrelevant to the target, and the 19 that matter are thin. The spline borrows strength across ages, so all 1,274 patients inform the age adjustment and all 27 exposed patients inform the target.

Understanding notes

Deviance. The handout now defines it inline, in the instructor’s 2026-10-04 edit, as −2 × log-likelihood, which is exact for 0/1 binomial data because the per-observation saturated log-likelihood is zero (862.4044 by hand = -2*logLik(fit)), and AIC = deviance + 2p for every binomial fit here (866.40; 824.60; 825.95; 906.23). Expect the question “why isn’t a saturated model’s deviance zero?”: the deviance’s reference model has one parameter per row (1,274), not per covariate pattern; the same exposure-only model fit to the two-row grouped table (cbind(d, n - d) ~ early_vasopressor) has deviance 3e-13, and the age-saturated model on the 73-row grouped table has deviance ~0 with null deviance 105.24 = 865.47 − 760.23. Poisson caveat: for 0/1 “counts” the saturated log-likelihood is −Σy = −136, so Step 2’s Poisson deviances (863.03 offset, 592.88 free) are −2ℓ − 272 and AIC ≠ deviance + 2p there; the handout only compares deviances within one family on the same rows, so nothing bites.
The delta-method aside (added 2026-10-04 at the instructor’s request, as a collapsed optional box; simplified the same day to: one SE re-expressed → delta method defined → logistic vcov to log-RR SE to RD SE → other directions left for later; the hand formulas below are now key-only). First appearance of the delta method in the course (not in the syllabus or any lecture week). Gradients: with $p = \text{expit}(\eta)$, $dp/d\eta = p(1-p)$ and $d\log p/d\eta = 1 - p$; the exposed have $\eta_1 = \beta_0 + \beta_1$, the unexposed $\eta_0 = \beta_0$, so $\partial(p_1 - p_0)/\partial\beta = (p_1(1-p_1) - p_0(1-p_0),\; p_1(1-p_1))$ and $\partial(\log p_1 - \log p_0)/\partial\beta = (p_0 - p_1,\; 1 - p_1)$; then $\text{Var} = g' V g$. Hand formulas: binomial variances for the RD, $(1-p)/x$ per group for the log RR, Woolf’s $1/a+1/b+1/c+1/d$ for the log OR. Because the model is saturated, the Fisher information under any link is the information for $(p_0, p_1)$ transformed by that link’s Jacobian, so the three glm() variances are exact re-expressions of one another; the delta method recovers them to ~3e-6 relative, the residual being glm()’s last-iteration working weights (Woolf 0.472095 vs. glm 0.472093). The reverse direction works too (identity fit → log OR SE 0.4721). msm::deltamethod is not installed in the course environment, hence $g' V g$ by hand; a finite-difference gradient gives the same numbers. In an unsaturated model the links no longer share a likelihood maximum, so there is no exact correspondence to recover.
Framing (instructor, 2026-10-04). The link is explained through Greenland’s lossy-compression image, not as an “assumption”: a saturated model is lossless, so the link only picks the coordinates; an unsaturated model is lossy, so the link decides which pattern survives. Step 3 then varies how much is kept. Keep that language if editing; the quotations are from Greenland, Stat Med 2025;44:e10313, section 1.1.
No randomness. Every number is a deterministic function of the canonical Lab 1 dataset. If anything drifts, the dataset changed.
Step 1a converges under all three links without starting values because the saturated MLE (0.104, 0.222) is interior to the parameter space. The exposure-only identity fit is the one binomial identity model in the course that “just works”; students who add age to it reproduce the Step 1b failure.
Step 1b identity failure is a boundary, not a bug. Zero deaths at ages 18–27 (52 patients) put the MLE at risk = 0 for age 18. With start = c(0.001, 0.0017) and maxit = 500 glm reports convergence after 45 iterations at (−0.04122, 0.00229), fitted risk 0.000000 at age 18; different platforms may stop at a slightly different boundary point, which is why the handout draws the identity line from lm() instead and only quotes the error message.
tryCatch, not try. The handout captures the error message as a returned string so it prints in the rendered HTML and the chunk continues; try() would send the message to stderr, which knitr may or may not show.
The two-row SE quirk (Step 2). Coefficients agree to 2.6e-8, but summary.glm computes standard errors from the working weights of the last IRLS iteration, and the 1,274-row and 2-row fits take slightly different paths, so the exposure SE prints as 0.417518 vs. 0.417563. The handout therefore claims equality of coefficients only. (The no-offset individual fit’s SE, 0.4175631, matches the aggregated one exactly.)
The freed-coefficient deviance (592.9) is dramatically lower than the offset model’s (863.0): log follow-up is a near-outcome, which is the whole Q5 point. Early drafts of this handout quoted a smaller gap; the verified numbers are 863.03 / 592.88.
Predicting a factor(age) fit works only at observed ages (the figure uses sat_grd); a continuous grid errors with “factor has new levels”. glm stores xlevels, so a subset of observed ages is fine.
Tibble rounding. summarise() output prints with 3 significant digits (7310.06 shows as 7310., 820.60 as 821.), so the handout wraps the two tables students read numbers from in as.data.frame() / data.frame(). If a student’s table “disagrees” at the third digit, this is usually why.
Age 90 is a pooled stratum (everyone over 89, n = 59, Lab 1’s de-identification recode), so the saturated fit’s last point rests on a decade of ages and is its most stable “single-year” estimate; the spline’s downturn at 90 (0.201 vs. the linear fit’s 0.218) is the natural-spline boundary constraint meeting that pooled point.
The saturated fit emits no glm warning (fitted values of 8.6e-9 do not trip the eps test), so the only visible symptom of the 30 empty strata is the coefficient table: intercept −18.57 with SE 2,465. If a future R version changes the convergence tolerance, the intercept and SE will shift while the fitted values stay at 0.
Exposure OR under age-only adjustment is 2.16, not Lab 3’s 1.54. The handout says so explicitly (adjustment set is deliberately age alone). Expect the question anyway; the answer is sex, MAP, and HR.