Lab 5 — Fit the GLM Family to a Single Outcome

EPID 785R · Regression and Study Design

Author

YOUR NAME HERE

Published

October 19, 2026

Introduction

The last lecture made a claim about generalized linear models: there is one fitting loop, and every member of the family is that loop with two functions plugged in, a variance function and a link. The lecture closed with an experiment, one mean model fit under several distributions, which returned the same coefficient when the model was saturated and nearly the same coefficient when it was not.

In this lab, we explore that idea more in practice. We take the model from Lab 3 (in-hospital death against early vasopressor initiation, age, sex, baseline MAP, and baseline heart rate) and fit it as a logistic model, a log-binomial model, a Poisson model with and without the person-time offset, and a linear probability model: one left-hand side, one right-hand side, four distribution-and-link pairs. We watch which fits converge and which do not, what each one reports for the exposure, and where each one’s fitted values leave the probability scale. We then turn one more dial, the flexibility of the right-hand side, in an exercise carried over from the Lab 4 handout (which we did not get to; the handout is on the course site). Finally we open the fitted logistic model and look at what the last iteration of the fitting loop leaves behind: residuals, leverage, and Cook’s distance.

In this lab, we will:

  1. fit four models that share one right-hand side (logistic, log-binomial, Poisson with and without the person-time offset, and linear probability), noting for each its distribution, its link, and the contrast its exposure coefficient reports, and record what the fitting loop did: how many iterations, whether it converged, and what it needed to get there (nothing, starting values, or more iterations);
  2. compare their predicted probabilities, and find where each model’s fitted values leave the probability scale;
  3. fit one relationship (death against age) three ways: linear, natural cubic spline, and saturated (the nonparametric maximum likelihood estimator); visualize the spectrum they span; and watch what the choice does to the exposure coefficient sitting next to the age term;
  4. compute the logistic model’s residuals, leverage, and Cook’s distance, first by hand from the last iteration of the loop and then with R’s helper functions.

Deliverables (yours to keep, nothing is submitted): this document rendered to HTML with your answers in the answer blocks. Six 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 6 notes (Generalized Linear Models) and Clark & Berry’s Models Demystified, chapter 8 (reading list). Zou (2004) is behind the Poisson route to a risk ratio in Step 1; Greenland’s (2025) image/metaphor of a fitted model as a lossy compression of the data organizes Step 3; Pregibon (1981) is the source of the diagnostics in Step 4.

1 Before you start (~3 min)

  1. Open lab5.Rproj in RStudio.
  2. The lab needs your Lab 1 analytic dataset. As in Labs 2 to 4, a copy ships in this lab’s data/ folder in the starter package.
  3. Packages: the tidyverse, splines (ships with R; the natural cubic spline in Step 3), and sandwich + lmtest (Lab 2’s robust standard errors, used in two lines of Step 1).
pacman::p_load(tidyverse, splines, 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 the website into labs/lab5/data/.")
}
lab1 <- read_csv(path, show_col_types = FALSE)
nrow(lab1)   # expect 1,292

# As in Labs 2 to 4: 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

# Lab 3's fit used 1,271 of these rows: glm() silently dropped the three
# patients with no baseline heart rate. Form that sample explicitly, so
# that every model in this lab is fit to exactly the same rows.
dat <- coh |> drop_na(baseline_hr)
nrow(dat)    # expect 1,271

2 Step 1 — Fit the four (~16 min)

2.1 The logistic model

Lab 3’s model, refit (provided). The formula object f holds the right-hand side that every fit in this lab shares. Besides the printout, ask the fitted object what the loop did: how many iterations, and whether it converged.

# the one right-hand side, shared by every model in this lab
f <- hospital_death ~ early_vasopressor + age + sex + baseline_map + baseline_hr

fit_lgt <- glm(f, family = binomial(link = "logit"), data = dat)
round(summary(fit_lgt)$coefficients, 4)

c(iterations = fit_lgt$iter, converged = fit_lgt$converged)
round(exp(coef(fit_lgt)[["early_vasopressor"]]), 3)   # the odds ratio
round(range(fitted(fit_lgt)), 4)                      # every fitted risk

Checkpoints: Lab 3’s printout: the exposure coefficient is 0.429 (SE 0.509), odds ratio 1.54. The loop needed 6 iterations and converged. The fitted risks run from 0.005 to 0.691, inside (0, 1) as the logit link guarantees.

2.3 The Poisson route to the risk ratio

The lecture’s first routine repair: keep the log link, swap the binomial for the Poisson, and replace the model-based standard error with a robust one (Zou 2004; Lab 2 did exactly this for the full-cohort risk ratio). Complete TODO (1b): the family.

# TODO (1b): the Poisson distribution with its canonical link
fit_pois <- glm(f, family = ???, data = dat)
c(iterations = fit_pois$iter, converged = fit_pois$converged)
round(exp(coef(fit_pois)[["early_vasopressor"]]), 3)   # the risk ratio
c(max_fitted = max(fitted(fit_pois)), above_1 = sum(fitted(fit_pois) > 1))

# the robust standard error, as in Lab 2 (tomorrow's lecture explains why)
round(coeftest(fit_pois, vcov = vcovHC(fit_pois, type = "HC3"))["early_vasopressor", ], 4)
round(sqrt(vcov(fit_pois)["early_vasopressor", "early_vasopressor"]), 4)   # model-based, for comparison

Checkpoints: 6 iterations, no starting values needed. Risk ratio 1.41 (the log-binomial said 1.29). One patient (the same one) gets a fitted “risk” of 1.006: the Poisson likelihood is defined there, so the loop never noticed. The exposure coefficient’s standard error is 0.425 model-based and 0.405 robust.

2.4 The rate model: the same Poisson machine with an offset

Add the log of each patient’s follow-up as an offset, a term whose coefficient is fixed at one, and the same Poisson model reports a rate ratio. Complete TODO (1c).

# TODO (1c): the log of person-time as an offset
fit_rate <- glm(hospital_death ~ early_vasopressor + age + sex + baseline_map +
                  baseline_hr + offset(???),
                family = poisson(link = "log"), data = dat)
c(iterations = fit_rate$iter, converged = fit_rate$converged)
round(exp(coef(fit_rate)[["early_vasopressor"]]), 3)   # the rate ratio

# what this model's fitted values are
c(min = min(fitted(fit_rate)), max = max(fitted(fit_rate)),
  above_1 = sum(fitted(fit_rate) > 1), total = sum(fitted(fit_rate)))

Checkpoints: 7 iterations; rate ratio 1.38. The fitted values are not risks: each is the number of deaths the model expects over that patient’s own follow-up, from 0.0001 to 2.95, seven of them above one, and they sum to 136.0, exactly the number of deaths. (That last equality is the canonical-link property from the lecture: with the log link and the Poisson distribution, the observed sum equals the fitted sum.)

2.5 The linear probability model

Last, the Gaussian distribution with the identity link, which is the ordinary linear model, and whose exposure coefficient is a risk difference. Complete TODO (1d): the family.

# TODO (1d): the Gaussian distribution with its canonical link
fit_lpm <- glm(f, family = ???, data = dat)
c(iterations = fit_lpm$iter, converged = fit_lpm$converged)
round(coef(fit_lpm)[["early_vasopressor"]], 4)     # the risk difference
max(abs(coef(fit_lpm) - coef(lm(f, data = dat))))   # the same as lm()?
c(min = min(fitted(fit_lpm)), below_0 = sum(fitted(fit_lpm) < 0),
  above_1 = sum(fitted(fit_lpm) > 1))

# standard errors: model-based, and robust (as in Lab 2)
round(c(model  = sqrt(vcov(fit_lpm)["early_vasopressor", "early_vasopressor"]),
        robust = sqrt(vcovHC(fit_lpm, type = "HC3")["early_vasopressor",
                                                     "early_vasopressor"])), 4)

# the binomial distribution with the identity link (provided): run it, read it
tryCatch(glm(f, family = binomial(link = "identity"), data = dat),
         error = function(e) conditionMessage(e))
tryCatch(glm(f, family = binomial(link = "identity"), data = dat,
             start = coef(fit_lpm)),
         error = function(e) conditionMessage(e))

Checkpoints: 2 iterations (the Gaussian-identity pair has a closed form; the loop’s second pass only confirms that nothing moved), and the coefficients equal lm()’s to machine precision. Risk difference 0.071 (model-based SE 0.059, robust 0.082). The smallest fitted risk is −0.13, and 106 patients get a negative one. The binomial-identity version refuses to start, even when handed the least-squares coefficients as a starting point: “cannot find valid starting values.” (The Lab 4 handout showed the one-covariate version of this refusal.)

2.6 Side by side

One row per fit (provided).

row_of <- function(fit) {
  b  <- coef(fit)[["early_vasopressor"]]
  se <- sqrt(vcov(fit)["early_vasopressor", "early_vasopressor"])
  data.frame(family = fit$family$family, link = fit$family$link,
             iterations = fit$iter, converged = fit$converged,
             estimate = round(b, 4), exp_estimate = round(exp(b), 3), SE = round(se, 4))
}
rbind(logistic       = row_of(fit_lgt),
      log_binomial   = row_of(fit_lb),
      poisson        = row_of(fit_pois),
      poisson_offset = row_of(fit_rate),
      linear_prob    = row_of(fit_lpm))

Checkpoints: exponentiated exposure coefficients 1.535, 1.285, 1.408, 1.384, and (not exponentiated) the risk difference 0.0713; iterations 6 / 29 / 6 / 7 / 2. (The exp_estimate column is meaningless for the last row.)

NoteQuestions — write your answers below each prompt

Q1. The log-binomial model failed at the default start, needed starting values and 29 iterations, and converged with one fitted risk at 0.997; without baseline heart rate it converged unaided in 6. Explain, in terms of the range of the fitted mean under each link, why a covariate running from 29 to 182 breaks this link and not the logit. Then: the lecture listed the routine repairs in order. Which risk ratio would we report, and what do we give up by reporting it?

YOUR ANSWER:

Q2. One exposure, four numbers: OR 1.54, RR 1.29 (log-binomial) or 1.41 (Poisson), rate ratio 1.38, RD 0.071. The two risk ratios come from models with the same link and the same right-hand side; only the distribution differs. In the lecture’s experiment, the crude risk ratio was identical across distributions and the adjusted one was not. Using what the distribution does inside the fitting loop, explain why 1.29 and 1.41 differ, and what would make them agree.

YOUR ANSWER:

Q3. The Poisson model’s model-based standard error (0.425) is about 5% larger than its robust one (0.405). In the lecture’s NHEFS example the model-based one was 60% too large. Why so different here? (Hint: compare the two variance functions, \(\mu\) and \(\mu(1-\mu)\), at \(\mu = 0.1\) and at \(\mu = 0.5\).)

YOUR ANSWER:

3 Step 2 — Predicted probabilities (~10 min)

Each of the four risk models assigns every patient a fitted value. Put them side by side (provided).

preds <- tibble(logistic     = fitted(fit_lgt),
                log_binomial = fitted(fit_lb),
                poisson      = fitted(fit_pois),
                linear_prob  = fitted(fit_lpm))

round(sapply(preds, function(p) c(min = min(p), median = median(p), mean = mean(p),
                                  max = max(p), below_0 = sum(p < 0),
                                  above_1 = sum(p > 1))), 4)
mean(dat$hospital_death)        # the observed risk
round(colSums(preds), 2)        # the fitted values summed, vs. 136 deaths
round(cor(preds), 3)

# the patient with the highest logistic risk, under all four models
top <- which.max(preds$logistic)
dat[top, c("early_vasopressor", "age", "sex", "baseline_map", "baseline_hr",
           "hospital_death")]
round(preds[top, ], 3)

Checkpoints: the mean fitted value is 0.1070 (the observed risk) for three of the four models and 0.1066 for the log-binomial; the sums are 136.00, 135.52, 136.00, 136.00. The medians of the first three agree to three decimals (0.089); the linear model’s is 0.107. The maxima are 0.691, 0.997, 1.006, 0.426; the linear probability model is the only one with fitted values below zero (106 of them). All four are highly correlated (0.99 among the first three; 0.91 between the logistic and the linear). The top patient is the one from Step 1 (76, male, MAP 31, heart rate 182, died): 0.691 / 0.997 / 1.006 / 0.426.

The figure (provided) holds that patient’s covariates fixed and lets heart rate run across the observed range, under each model.

ref <- data.frame(early_vasopressor = 0, age = 76, sex = "Male", baseline_map = 31,
                  baseline_hr = seq(40, 185, by = 1))
curves <- bind_rows(
  tibble(hr = ref$baseline_hr, risk = predict(fit_lgt,  ref, type = "response"),
         Model = "logistic"),
  tibble(hr = ref$baseline_hr, risk = predict(fit_lb,   ref, type = "response"),
         Model = "log-binomial"),
  tibble(hr = ref$baseline_hr, risk = predict(fit_pois, ref, type = "response"),
         Model = "Poisson (log link)"),
  tibble(hr = ref$baseline_hr, risk = predict(fit_lpm,  ref, type = "response"),
         Model = "linear probability")) |>
  mutate(Model = factor(Model, c("logistic", "log-binomial", "Poisson (log link)",
                                 "linear probability")))

ggplot() +
  geom_hline(yintercept = c(0, 1), linetype = "dotted", color = "grey50") +
  geom_rug(data = dat, aes(x = baseline_hr), sides = "b", alpha = .15) +
  geom_line(data = curves, aes(x = hr, y = risk, color = Model), linewidth = .8) +
  scale_color_manual(values = c("#0072B2", "#D55E00", "#009E73", "#CC79A7")) +
  labs(x = "Baseline heart rate (beats/min)",
       y = "Fitted risk of in-hospital death\n(unexposed man, 76, MAP 31)",
       color = NULL) +
  theme_classic() + theme(legend.position = "bottom")

The four curves are four compressors tuned to four patterns: a sigmoid (logit), an exponential (log, twice, with two sets of weights), and a straight line (identity). Over the range where most patients sit (indicated via the rug plot) they nearly coincide; they deviate where the data thin out. The exponential curves cross one near the top of the observed range (the Poisson’s at a heart rate of 182, the log-binomial’s at 183); the sigmoid bends toward one and stops at 0.69.

NoteQuestions

Q4. Three of the four models return fitted values that sum exactly to the 136 observed deaths; the fourth sums to 135.5. Which three, and which property from the lecture explains it? Why does the log-binomial model not share it?

YOUR ANSWER:

Q5. For one patient profile the four models agree over most of the heart-rate range and disagree at the top: 0.69, 0.997, 1.006, 0.43 at 182 beats per minute. The lecture’s closing table had a “where it breaks” column. Match each curve to its entry. Then say which model we would use to report a predicted risk for a new patient like this one, and why the answer could differ from the model we would use to report the exposure contrast.

YOUR ANSWER:

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

Steps 1 and 2 turned the distribution-and-link dial with the right-hand side held fixed. This step holds the distribution and link fixed (binomial, logit) and changes the flexibility encoded in the right-hand side.

Recall the image from the Week 4 and 5 lectures: a fitted model is a lossy compression of the data. In Greenland’s words, a parametric model form is a filter that passes along “information that follows the form” and filters out “information which does not so conform.” The flexibility of the right-hand side is how “lossy” the compression becomes: how many parameters the model spends and therefore how much of the pattern between the outcome and a covariate it keeps. The Week 5 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.

Age lets us see the whole spectrum in one relationship: 73 distinct values, with between 2 and 59 patients at each. That is enough to fit a model with a parameter for every observed year of age. 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 below 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.
fit_lin <- glm(hospital_death ~ age,             family = binomial, data = dat)
fit_spl <- glm(hospital_death ~ ns(age, df = 4), family = binomial, data = dat)
fit_sat <- glm(hospital_death ~ factor(age),     family = binomial, data = dat)

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
c(iterations_saturated = fit_sat$iter)

Checkpoints: 2, 5, and 73 parameters; residual deviances 820.25, 815.67, 759.76; AICs 824.25, 825.67, 905.76. (The null deviance, for reference, is 864.79.) 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. The saturated fit needed 17 iterations, three times the logistic fits of Step 1; the printout below says why.

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

by_age <- dat |>
  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

# where the saturated model is "estimating" from very little
by_age |> filter(age %in% c(18, 28, 44, 85, 88, 90))
sum(by_age$deaths == 0)                      # ages with no deaths at all

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

Checkpoints: the difference is zero to numerical precision (about 1e−8: the loop stops when it is close enough). The saturated model is the stratum-specific proportions. The filtered rows show what those proportions rest on: 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 (where Lab 1 placed everyone over 89). 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, while age 28, with one death, gets a coefficient of +16.4 to climb back. This is the lecture’s second kind of convergence failure, separation: not a fitted mean leaving its range (Step 1’s log-binomial) but a coefficient leaving the real line, produced here not by an unusual covariate but by asking for one parameter per age from a cohort that cannot supply it.

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(dat$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")

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.

4.1 The target next to the nuisance

Bring the exposure back, with the rest of Step 1’s right-hand side. Under Lab 3’s causal reading, the exposure coefficient is the target and age is part of the 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 the age term beside it changes and nothing else does (provided). The first row is Step 1’s logistic fit.

or_row <- function(fit) {
  b  <- coef(fit)[["early_vasopressor"]]
  se <- 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(fit_lgt),
  age_spline    = or_row(glm(hospital_death ~ early_vasopressor + ns(age, df = 4) +
                               sex + baseline_map + baseline_hr,
                             family = binomial, data = dat)),
  age_saturated = or_row(glm(hospital_death ~ early_vasopressor + factor(age) +
                               sex + baseline_map + baseline_hr,
                             family = binomial, data = dat))), 3)

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

Checkpoints: OR 1.54 (SE 0.509) with age linear, 1.56 (SE 0.508) with the spline, 1.72 (SE 0.545) with age saturated, on 6, 9, and 77 parameters; only 19 of the 73 age strata contain an exposed patient.

NoteQuestions

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

YOUR ANSWER:

Q7. Going from a linear to a spline age term moved the exposure OR from 1.54 to 1.56 and left its SE at 0.51; going to a saturated age term moved it to 1.72 and the SE to 0.55. Which change is providing bias protection in the nuisance function, and which is resulting in a variance cost? Why does saturation hurt here in particular? (Hint: the last line of the chunk, and the 27 exposed patients.)

YOUR ANSWER:

5 Step 4 — What the last iteration leaves behind (~16 min)

The fitting loop ends with one weighted least-squares fit: a working response \(z\), a diagonal weight matrix \(W\) (for the logit link, \(w_i = \hat\mu_i(1 - \hat\mu_i)\)), and the design matrix \(X\). Three things fall out of that last fit for free, and they are the diagnostics of this step.

  • Residuals, in four currencies. The response residual is the raw gap, \(y_i - \hat\mu_i\). The Pearson residual scales it by the standard deviation the model expects at that mean, \((y_i - \hat\mu_i) / \sqrt{V(\hat\mu_i)}\), so that its squares sum to the Pearson statistic. The deviance residual is the signed square root of observation \(i\)’s contribution to the deviance, so that its squares sum to the deviance. The working residual is the residual of the final weighted regression, on the scale of \(z\). For a binary outcome every one of these takes only two values at a given \(\hat\mu_i\) (one for a death, one for a survival), which is why a plot of residuals against fitted values looks nothing like the one we know from linear regression.
  • Leverage, the diagonal of the hat matrix of that final regression, \(H = W^{1/2} X (X'WX)^{-1} X' W^{1/2}\). In linear regression leverage depends on \(X\) alone; here \(W\) depends on \(\hat\mu\), so leverage depends on the outcome too. The diagonal sums to the number of parameters, so the average leverage is \(p/n\).
  • Cook’s distance, a one-number summary of how far the whole coefficient vector would move if observation \(i\) were dropped, approximated without refitting as \(D_i = \dfrac{r_{P,i}^2}{p} \cdot \dfrac{h_i}{(1 - h_i)^2}\) (Pregibon 1981): a large residual and high leverage.

Everything below is on Step 1’s logistic model.

5.1 Residuals

Complete TODO (4a): the response and Pearson residuals by hand, then check them against residuals().

mu <- fitted(fit_lgt)

# TODO (4a): the raw gap, and the raw gap scaled by the binomial standard deviation
r_resp <- ??? - mu
r_pear <- r_resp / sqrt(???)
r_dev  <- residuals(fit_lgt, type = "deviance")

c(response = max(abs(r_resp - residuals(fit_lgt, type = "response"))),
  pearson  = max(abs(r_pear - residuals(fit_lgt, type = "pearson"))))   # both zero

round(rbind(response = range(r_resp), pearson = range(r_pear), deviance = range(r_dev),
            working  = range(residuals(fit_lgt, type = "working"))), 3)
c(sum_dev_sq = sum(r_dev^2), deviance = deviance(fit_lgt), sum_pearson_sq = sum(r_pear^2))

# the largest Pearson residual belongs to...
dat[which.max(r_pear), c("early_vasopressor", "age", "sex", "baseline_map",
                         "baseline_hr", "hospital_death")]
round(mu[[which.max(r_pear)]], 4)

Checkpoints: both by-hand residuals match R’s exactly. The response residuals run from −0.572 to 0.976, the Pearson from −1.156 to 6.323, the deviance from −1.303 to 2.725, and the working from −2.3 to 41. The squared deviance residuals sum to 781.62, the model’s deviance; the squared Pearson residuals sum to 1,249.4. The largest Pearson residual belongs to a 28-year-old man (MAP 62, heart rate 101) who died with a fitted risk of 0.024: the only death among the 146 patients under 40.

Two plots (provided). The first is the residual-against-fitted plot we know from linear regression, which carries no information here: a binary outcome gives it two bands by construction. The second, Gelman and Hill’s binned residual plot, is the one to read: patients are grouped into deciles of fitted risk, and each decile’s mean residual is drawn with a ±2 standard-error band for a well-calibrated model.

ggplot(data.frame(mu, r_dev, death = factor(dat$hospital_death)),
       aes(x = mu, y = r_dev, color = death)) +
  geom_point(alpha = .4, size = 1) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("survived", "died")) +
  labs(x = "Fitted risk", y = "Deviance residual", color = NULL) +
  theme_classic() + theme(legend.position = "bottom")
bins <- dat |>
  mutate(p = mu, r = r_resp, bin = ntile(p, 10)) |>
  group_by(bin) |>
  summarise(n = n(), p_mean = mean(p), observed = mean(hospital_death),
            r_mean = mean(r), se = sqrt(p_mean * (1 - p_mean) / n), .groups = "drop")
as.data.frame(round(bins, 4))

ggplot(bins, aes(x = p_mean, y = r_mean)) +
  geom_ribbon(aes(ymin = -2 * se, ymax = 2 * se), fill = "grey85") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_point(size = 2) + geom_line() +
  labs(x = "Mean fitted risk in decile", y = "Mean residual (observed − fitted)") +
  theme_classic()

Checkpoints: the first plot is two bands, deaths above zero and survivals below, each a smooth function of the fitted risk; it carries no information about fit beyond the two labels. In the binned plot nine of the ten deciles sit inside the band; the second decile (mean fitted 0.035, observed 0.071) sits just outside it.

5.2 Leverage

Complete TODO (4b): the hat diagonal by hand from the last iteration’s weights and the design matrix, checked against hatvalues().

X <- model.matrix(fit_lgt)
# TODO (4b): the working weights of the last iteration, stored in the fitted object
w <- ???
max(abs(w - mu * (1 - mu)))                                   # they are mu(1 - mu)

# diag of W^1/2 X (X'WX)^-1 X' W^1/2, one row at a time
h_hand <- rowSums(((sqrt(w) * X) %*% solve(t(X) %*% (w * X))) * (sqrt(w) * X))
h      <- hatvalues(fit_lgt)
max(abs(h_hand - h))                                          # zero
c(sum = sum(h), mean = mean(h), p_over_n = length(coef(fit_lgt)) / nrow(dat),
  max = max(h))

# who carries the most leverage?
top_h <- order(-h)[1:6]
cbind(dat[top_h, c("early_vasopressor", "age", "sex", "baseline_map", "baseline_hr",
                   "hospital_death")],
      leverage = round(h[top_h], 4), fitted = round(mu[top_h], 3))
round(tapply(h, dat$early_vasopressor, mean), 4)   # mean leverage, unexposed vs. exposed

Checkpoints: the stored weights are \(\hat\mu(1 - \hat\mu)\) (to 1e−8), the by-hand diagonal matches hatvalues() to machine precision, and it sums to 6, the number of parameters, for a mean of 0.0047. The maximum is 0.0765, sixteen times the mean. All six top-leverage patients are exposed, and each also sits at an edge of the design in a second direction: a MAP at the bottom of its range (32 to 38) or at the top (64), or an age of 89 or 90; the mean leverage of an exposed patient is 0.043, eleven times an unexposed patient’s 0.004.

5.3 Cook’s distance

Provided: R’s cooks.distance(), the by-hand formula, the six most influential patients, and the refit without the most influential one.

cd      <- cooks.distance(fit_lgt)
cd_hand <- (r_pear^2 / length(coef(fit_lgt))) * h / (1 - h)^2
max(abs(cd_hand - cd))                                        # zero

top_c <- order(-cd)[1:6]
cbind(dat[top_c, c("early_vasopressor", "age", "sex", "baseline_map", "baseline_hr",
                   "hospital_death")],
      leverage = round(h[top_c], 4), cooks = round(cd[top_c], 4),
      fitted = round(mu[top_c], 3))

# refit without the single most influential patient
fit_drop1 <- glm(f, family = binomial, data = dat[-top_c[1], ])
round(rbind(all_1271 = coef(fit_lgt), drop_one = coef(fit_drop1)), 4)
round(exp(c(all_1271 = coef(fit_lgt)[["early_vasopressor"]],
            drop_one = coef(fit_drop1)[["early_vasopressor"]])), 3)

# for tomorrow: the linear probability model's residual variance, by tertile of fitted risk
round(tapply(residuals(fit_lpm), ntile(fitted(fit_lpm), 3), var), 4)

Checkpoints: the by-hand Cook’s distances match R’s. The six most influential patients are five of the six exposed patients who died (fitted risks 0.16 to 0.28, Cook’s 0.034 down to 0.025) and the 28-year-old from the residual step (the sixth exposed death is seventh). Dropping the single most influential patient moves the exposure coefficient from 0.429 to 0.202 (OR 1.54 to 1.22) and every other coefficient by less than 0.06. The last line is for tomorrow’s lecture: the linear probability model’s residual variance grows from 0.041 to 0.066 to 0.160 across tertiles of fitted risk, in a model whose distribution says it is constant. That is the binomial variance function \(\mu(1 - \mu)\) showing through, and the robust standard error we computed in Step 1 is the repair.

NoteQuestions

Q8. The six highest-leverage patients are all exposed, and the mean leverage of an exposed patient is eleven times that of an unexposed one. Using the formula for \(H\) and the fact that only 27 of 1,271 patients are exposed, explain why. Then explain, in one or two sentences, why leverage in a GLM depends on the outcome when in linear regression it does not.

YOUR ANSWER:

Q9. Cook’s distance says the exposure coefficient rests on a handful of patients: remove one exposed death and the OR drops from 1.54 to 1.22. Is this a flaw in the model, a flaw in the data, or neither? What would we report about this analysis, and how does what Cook’s distance found relate to something the printout already told us in Lab 3 (SE 0.51, p = 0.40)?

YOUR ANSWER:

References

  • Clark M, Berry S. Models Demystified: A Practical Guide from Linear Regression to Deep Learning, chapter 8 (Generalized Linear Models). https://m-clark.github.io/book-of-models/
  • Zou G. A Modified Poisson Regression Approach to Prospective Studies with Binary Data. Am J Epidemiol. 2004;159(7):702–706. (The Poisson route to the risk ratio, Step 1.)
  • Marschner IC, Gillett AC. Relative Risk Regression: Reliable and Flexible Methods for Log-Binomial Models. Biostatistics. 2012;13(1):179–192. (Dedicated algorithms for the log-binomial fit.)
  • Greenland S. Some Ways to Make Regression Modeling More Helpful Than Misleading. Stat Med. 2025;44:e10313. (The lossy-compression image, Step 3; in the course PDF of Carlin & Moreno-Betancur.)
  • Greenland S, Mansournia MA, Altman DG. Sparse Data Bias: A Problem Hiding in Plain Sight. BMJ. 2016;352:i1981. (Background for Q7.)
  • Pregibon D. Logistic Regression Diagnostics. Ann Stat. 1981;9(4):705–724. (Leverage and Cook’s distance for GLMs, Step 4.)
  • Gelman A, Hill J. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, 2007, section 5.6. (The binned residual plot.)
  • Fox J. Applied Regression Analysis and Generalized Linear Models, 3rd ed. Sage, 2016, chapter 15. (GLM diagnostics.)
  • Pollard T, et al. eICU Collaborative Research Database Demo (v2.0.1). PhysioNet, 2021. doi:10.13026/4mxk-na84 (the Lab 1 source data).