Lab 5 — Instructor answer key

Companion file: solution.R (in this folder), the lab chunk code with the six TODOs (1a, 1b, 1c, 1d, 4a, 4b) completed, 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 (the alternative start vector, HC0 versus HC3, the rate model’s robust SE, the decade-bin model, fitted risks at age 61, exact intervals for ages 85 and 88, the refit without all six influential patients, and the separation that follows from dropping every exposed death).

Contents: TODO code answers · the four-model table (Step 0) · expected checkpoints · model answers Q1–Q9 · understanding notes

TODO code answers

TODO Answer
1a st <- c(log(mean(dat$hospital_death)), rep(0, 5)), i.e. the intercept at log(136/1271) = −2.2349 and the five slopes at zero, so that every patient’s starting fitted risk is the marginal risk 0.107, inside (0, 1). Students who write log(0.107) by hand get −2.235 and the same fit. With the default maxit = 25 the loop stops unconverged ($converged is FALSE); with maxit = 500 it converges at iteration 29. Both calls print step size truncated warnings: the loop is halving its steps to stay inside the parameter space.
1b family = poisson(link = "log"), or simply family = poisson (log is the Poisson’s canonical link, so it is the default). Converges in 6 iterations with no starting values; one fitted value exceeds 1 (1.006) because the Poisson likelihood is defined there.
1c offset(log(followup_days)). Equivalent spellings: the offset = log(followup_days) argument to glm(). Students who write offset(followup_days) (no log) get a different model (an additive offset on the log scale) with a rate ratio of 1.40 and should be sent back to the Week 5 definition.
1d family = gaussian(link = "identity"), or family = gaussian. Two iterations; identical to lm(f, data = dat) to machine precision. Students who try binomial(link = "identity") here meet the error the provided lines then show on purpose.
4a r_resp <- dat$hospital_death - mu and r_pear <- r_resp / sqrt(mu * (1 - mu)). Both match residuals(fit_lgt, type = ...) exactly (max difference 0).
4b w <- fit_lgt$weights, the working weights of the last IRLS iteration. They equal mu * (1 - mu) to 1.8e−8 (the last iteration’s weights are computed from the previous iteration’s coefficients). Students who write w <- mu * (1 - mu) get a hat diagonal that matches hatvalues() to ~1e−9 instead of 1e−15, which is fine.

The four models’ parts (the former Step 0)

The handout’s Step 0, a parts-list table for students to complete, was cut by the instructor on 2026-10-13 for time; the formula f it defined now sits at the top of Step 1’s logistic chunk. The table is kept here because Q1 and Q4 draw on it. Model answers for the table; any phrasing that names the right object is fine. “Canonical” means the link under which the linear predictor is the distribution’s natural parameter (lecture, section 3).

model distribution link canonical? the exposure coefficient is a …
logistic binomial logit yes log odds ratio; exp() = the odds ratio of death comparing early initiators to non-initiators of the same age, sex, MAP, and heart rate
log-binomial binomial log no log risk ratio; exp() = the conditional risk ratio
Poisson, no offset Poisson log yes log risk ratio again: with a 0/1 outcome and no offset the Poisson mean is the risk, so exp() = the conditional risk ratio (Zou 2004); the distribution is “wrong” (variance μ, not μ(1 − μ)), which is why the standard error is replaced by a robust one
Poisson with offset(log(followup_days)) Poisson log yes log rate ratio; exp() = the conditional ratio of deaths per person-day; the offset’s coefficient is fixed at one
linear probability Gaussian identity yes risk difference (not exponentiated): the conditional difference in the probability of death; the distribution is “wrong” (constant variance for a 0/1 outcome), again repaired with a robust SE

Expected checkpoints

Verified 2026-10-13 and re-verified 2026-10-19 against the canonical Lab 1 dataset (1,292 rows → 1,274 with an outcome → 1,271 with a baseline heart rate; every model uses the 1,271) on R 4.6.1. No seeds, no randomness anywhere.

setup:   nrow(lab1) = 1292; nrow(coh) = 1274; nrow(dat) = 1271 (136 deaths, 27 exposed)
f + links: binomial()$link logit; poisson()$link log; gaussian()$link identity
Step 1:  logistic     iter 6, converged; exposure 0.4286 (SE 0.5090), OR 1.535; fitted 0.0048-0.6906;
                      deviance 781.62, AIC 793.62 (Lab 3's numbers)
         log-binomial default start -> "no valid set of coefficients has been found: please supply starting values"
                      start = (-2.2349, 0,0,0,0,0): maxit 25 -> iter 25, converged FALSE
                        (warnings: step size truncated due to divergence / out of bounds; algorithm did not converge)
                      maxit 500 -> iter 29, converged TRUE; exposure 0.2509 (SE 0.3635), RR 1.285;
                        max fitted 0.9974 (one patient above 0.99); deviance 782.12, AIC 794.12
                      key-only: start = (-5, 0...) -> iter 34, converged, RR 1.2852 (same answer)
                      without baseline_hr: iter 6, converged, RR 1.491, max fitted 0.422
                      boundary patient: unexposed, 76, Male, MAP 31, HR 182, died; logistic 0.691 / log-binomial 0.997
                      heart-rate range 29-182
         Poisson      iter 6; exposure 0.3419, RR 1.408; max fitted 1.0059 (1 above 1)
                      SE model 0.4254 / robust HC3 0.405 (HC0 0.385) / logistic's 0.509
         rate model   iter 7; exposure 0.3250, IRR 1.384; fitted 0.0001-2.9498, 7 above 1, sum 136.00
                      key-only: robust HC3 SE 0.5319 vs model 0.4256
         LPM          iter 2; RD 0.0713; = lm() exactly; fitted -0.1316..0.4263; 106 below 0, 0 above 1
                      SE model 0.0590 / robust HC3 0.0823
                      binomial identity: "no valid set of coefficients..."; with start = coef(fit_lpm):
                        "cannot find valid starting values: please specify some"
         side by side: exp(estimate) 1.535 / 1.285 / 1.408 / 1.384 / (1.074, meaningless; RD 0.0713)
                      iterations 6 / 29 / 6 / 7 / 2; SE 0.5090 / 0.3635 / 0.4254 / 0.4256 / 0.0590
Step 2:  means 0.1070 / 0.1066 / 0.1070 / 0.1070 (observed 0.1070); medians 0.0890 / 0.0889 / 0.0887 / 0.1069
         sums 136.00 / 135.52 / 136.00 / 136.00; maxima 0.6906 / 0.9974 / 1.0059 / 0.4263; LPM below 0: 106
         correlations with the logistic: 0.987 / 0.987 / 0.909 (log-binomial vs Poisson 0.999)
         top patient (row 526 of dat): 0.691 / 0.997 / 1.006 / 0.426
         key-only: lowest LPM fitted -0.1316 (19-year-old woman, MAP 60, HR 60; logistic 0.0048);
                   negative LPM risks reach age 69; all 106 are unexposed
         figure (unexposed man, 76, MAP 31): at HR 60 0.111/0.105/0.106/0.130; at 120 0.340/0.318/0.321/0.276;
                   at 182 0.691/0.997/1.006/0.426; log-binomial crosses 1 at HR 183, Poisson at 182
Step 3:  linear 2 params, deviance 820.25, AIC 824.25
         ns(age, 4) 5 params, deviance 815.67, AIC 825.67 (knots 55/68/78, boundary 18/90)
         factor(age) 73 params, deviance 759.76, AIC 905.76 (17 iterations, converged)
         null deviance 864.79
         LRT spline vs linear 4.59 on 3 df, p = 0.205; saturated vs spline 55.9 on 68 df, p = 0.85
         max |p_obs - saturated fit| 8.6e-9
         by age: 18 -> 0/7; 28 -> 1/10; 44 -> 0/3; 85 -> 6/15; 88 -> 1/13; 90 -> 14/59
         30 ages (253 patients) have zero deaths; n per age min 2, median 15, max 59
         saturated printout: intercept -18.57 (SE 2465); factor(age)19 0.00 (SE 4088); factor(age)28 +16.37
         key-only: decade bins (<30, 30s, ..., 80+) 7 params, deviance 814.21, AIC 828.21
                   fitted at 61: linear 0.079, spline 0.080, saturated 0.077 (1 of 13); saturated at 60 = 0 (0 of 22)
                   exact 95% CI: age 85 (6/15) 0.16-0.68; age 88 (1/13) 0.00-0.36
                   spline fitted at 18/40/70/90: 0.001/0.034/0.091/0.201; linear 0.015/0.035/0.110/0.218
         target next to nuisance (Step 1's RHS): OR 1.535 (SE 0.509, CI 0.57-4.16, 6 params)
                                                  1.558 (SE 0.508, CI 0.58-4.22, 9 params)
                                                  1.724 (SE 0.545, CI 0.59-5.02, 77 params)
         19 of 73 age strata contain an exposed patient
Step 4:  residuals by hand = residuals() exactly (0 difference)
         ranges: response -0.572..0.976; Pearson -1.156..6.323; deviance -1.303..2.725; working -2.337..40.978
         sum dev^2 781.62 = deviance; sum Pearson^2 1249.43
         largest Pearson residual: row 969, unexposed man, 28, MAP 62, HR 101, died, fitted 0.0244
         binned (deciles): mean residuals -0.008, +0.036, -0.010, -0.017, 0.000, -0.002, -0.021, -0.015, +0.038, +0.001;
                   SEs 0.011-0.041; one decile (the 2nd: fitted 0.035, observed 0.071) outside +/-2 SE
         key-only, by age band: 30-39 observed 0 of 73 vs fitted 0.028 (1.4 SE); all other bands inside
         leverage: fit$weights = mu(1-mu) to 1.8e-8; by hand = hatvalues() to 2.7e-15
                   sum 6, mean 0.00472 = 6/1271, max 0.0765; 95 patients above 2p/n, 55 above 3p/n
                   top six (all exposed): leverage 0.0765 / 0.0712 / 0.0706 / 0.0681 / 0.0625 / 0.0608,
                     profiles M81 MAP32 / M78 MAP64 / M67 MAP36 / F78 MAP38 / M89 MAP53 (died) / M90 MAP59
                   mean leverage unexposed 0.0039 vs exposed 0.0426; 19 of the 27 exposed sit in the top 27
         Cook's: by hand = cooks.distance() to 6.9e-18; max 0.0339; 121 patients above 4/n
                   top six: F72 (0.0339), F78 (0.0313), M66 (0.0288), M67 (0.0256), F78 (0.0251), all exposed deaths
                     with fitted 0.156-0.283; then the 28-year-old (0.0179, unexposed, fitted 0.024);
                     the six exposed deaths' ranks are 1 2 3 4 5 7
                   drop the top patient: exposure 0.4286 -> 0.2017, OR 1.535 -> 1.223; other coefs move <= 0.058
                   key-only: drop all six -> OR 0.213; drop the six exposed deaths -> coef -14.89, SE 487 (separation)
         LPM residual variance by tertile of fitted risk: 0.0411 / 0.0664 / 0.1600

Model answers, Q1–Q9

Numbered as in the handout after the instructor’s 2026-10-13 trim, which cut the original Q7.

Q1. Why heart rate breaks the log link and not the logit; what to report. Under the logit link the fitted risk is expit(η), which lies in (0, 1) for every η, so no covariate value can push a fitted mean out of the binomial’s range: the loop always has a likelihood to evaluate. Under the log link the fitted risk is exp(η), bounded below by zero but not above by one, and η rises by 0.018 per beat per minute. A patient at the top of a covariate with a 150-unit range (29 to 182), combined with age 76 and a MAP of 31, has η at or above zero, i.e. a fitted risk at or above one, where the binomial log-likelihood (log(1 − μ) for a survivor) does not exist. From the default start the loop cannot take a first step; from the marginal-risk start it walks toward that boundary and has to halve its steps to stay inside (“step size truncated”), which is why it needs 29 iterations, and it stops with one fitted risk at 0.997. Without heart rate the linear predictor never gets near zero (max fitted 0.42), so the constraint never binds and the loop converges unaided in 6. Which to report: the lecture’s order was Poisson-log with a robust variance (good) or standardization from the logistic fit (better, Week 8). The Poisson route gives RR 1.41 with robust SE 0.405 and is the honest default here: it converges without help and its standard error is right once made robust. What is given up: the binomial likelihood (one patient’s fitted “risk” is 1.006, and the estimate differs from the log-binomial’s 1.29, so the number depends on the route), and a little efficiency. Reporting the log-binomial’s 1.29 means reporting a fit coaxed to the boundary, whose existence depends on the start and on the sample, which the handout’s “reassuring, not a guarantee” sentence is about. Reporting the logistic’s OR 1.54 means reporting a different contrast, which is only a problem if a risk ratio was the question.

Q2. Why 1.29 and 1.41 differ when only the distribution differs. In the fitting loop the distribution enters only through the variance function in the weights, \(w_i = (d\mu_i/d\eta_i)^2 / V(\mu_i)\). With a log link \(d\mu/d\eta = \mu\), so the binomial weight is \(\mu^2/[\mu(1-\mu)] = \mu/(1-\mu)\) and the Poisson weight is \(\mu^2/\mu = \mu\). Same link, same working response, different weights: the two fits are two differently weighted compressions of the same data, and they differ wherever the weights disagree, which is wherever \(\mu\) is not small, exactly the high-risk patients who decide the fit at the top of the heart-rate range. When the model is saturated (the lecture’s crude experiment; Lab 2’s exposure-only Poisson, RR 2.13 under either distribution) the weights do not matter: the fitted values are the observed proportions whatever the weights, nothing is lost, and the two routes agree exactly. Here the compression is lossy, so the weights choose what is kept, and the two answers differ. The gap is larger than the lecture’s (1.33 vs 1.36 on NHEFS) because the log-binomial sits at its boundary: the constraint \(\mu \le 1\) binds for the HR-182 patient and bends the whole fit, whereas the Poisson is free to give him 1.006. What would make them agree: a saturated right-hand side; or a cohort in which every fitted risk is small, so that \(\mu/(1-\mu) \approx \mu\) for everyone (not this one: the top patients are at 0.5 to 0.7).

Q3. Why the Poisson SE is only 5% too large here but 60% on NHEFS. The Poisson variance function is \(\mu\) and the binomial’s is \(\mu(1-\mu)\); the ratio is \(1-\mu\). At \(\mu = 0.1\) the Poisson overstates the variance by a factor of 1/0.9 (11%, so the SE by about 5%); at \(\mu = 0.5\) by a factor of 2 (SE by 41%). The NHEFS outcome was above-median weight gain, a 50% outcome, so the Poisson variance function was badly wrong and the model-based SE 0.0795 was far above the robust 0.0520. Here the risk is 0.107 on average and below 0.15 for three-quarters of the cohort, so \(\mu \approx \mu(1-\mu)\) for almost everyone and the robust SE (0.405) barely corrects the model-based one (0.425). The repair is the same in both cases; how much it matters depends on how far the outcome’s variance is from the distribution’s variance function over the range of fitted means. (The log-binomial’s own SE, 0.364, is smaller still, since it uses the right variance function and a different, boundary-constrained fit.)

Q4. Which fitted values sum to 136, and why. The logistic (binomial-logit), the Poisson (Poisson-log), and the linear probability model (Gaussian-identity): the three canonical pairs. For a canonical link the score equations reduce to \(X'y = X'\hat\mu\), the observed sums equal the fitted sums for every column of \(X\), and the intercept column of ones makes that \(\sum y_i = \sum \hat\mu_i = 136\) (and, through the exposure column, 6 fitted deaths among the 27 exposed). The log link is not canonical for the binomial, so its score equations are \(\sum_i x_i (y_i - \hat\mu_i)/(1 - \hat\mu_i) = 0\), a weighted sum of residuals, and the plain sum need not be zero: 135.52. (Lab 4’s handout had the same fact from the other side: in the saturated exposure-only model every link reproduces the two observed proportions, because a saturated model has nothing to lose.)

Q5. Matching the curves to “where it breaks”; which model for a predicted risk. Logistic: its entry is separation, a coefficient leaving the real line, not a mean leaving its range; its curve bends toward one and stops at 0.69, and nothing in this cohort separated it (Step 3’s saturated age term is where separation appears). Log-binomial: fitted risk above 1; its curve reaches 0.997 at 182 and would cross one at 183, which is why the fit needed starting values and stopped at the boundary. Poisson: the lecture’s table says overdispersion, which is the count-outcome failure; for a binary outcome its breakage is the same range problem, met silently (1.006, no error, no warning) plus the wrong variance function (repaired by the robust SE). Linear probability: heteroscedasticity (Step 4’s last line) and fitted values outside [0, 1] (106 negatives; 0.43 here only because the line is too flat to reach the top). For a predicted risk for a new patient, the logistic model: its fitted values are probabilities everywhere, including at covariate values far from the center where an individual prediction is most wanted; it is calibrated in deciles (Step 4); and its likelihood has one maximum. The other three produce “risks” outside [0, 1] exactly in that region. For the exposure contrast the question chooses the scale (a risk difference, a risk ratio, an odds ratio), and the model that reports it directly with a valid standard error (the LPM with a robust SE for the RD; the Poisson with a robust SE for the RR) is a defensible choice even though its fitted values are not all probabilities, since the coefficient is the deliverable, not the fitted values; Week 8’s standardization from the logistic fit is the way to have both.

Q6. 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.6 of deviance (not quite the 6 they cost), and the saturated fit’s 68 extra parameters bought 55.9 (far short of 136). The likelihood-ratio 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 the target-next-to-nuisance table shows it does so at no cost in the target’s precision (SE 0.509 → 0.508), 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.

Q7. Bias protection vs. variance cost; why saturation hurts here. Linear → spline (1.54 → 1.56, SE unchanged at 0.51): the nuisance function buying bias protection, which in these data turned out to be nearly unnecessary (the age relationship is close to log-linear), at no cost in precision. Spline → saturated (1.56 → 1.72, SE 0.51 → 0.55): 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 patient 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,271 patients inform the age adjustment and all 27 exposed patients inform the target.

Q8. Why the exposed carry the leverage; why leverage depends on the outcome. The hat diagonal is \(h_i = w_i\, x_i'(X'WX)^{-1}x_i\): it is large when \(x_i\) is far from the weighted center of the design in a direction where the design is thin, because \((X'WX)^{-1}\) is large in that direction. The exposure column has 27 ones among 1,271 zeros, so every exposed row is far from the bulk in that direction and the matrix inverse is large there (the exposure’s SE, 0.51, is the largest in the table for the same reason). The six top rows combine exposure with a MAP at an end of its range (32 to 38, or 64) or an extreme age (89 to 90), far from the center in two directions at once. Hence a mean leverage of 0.043 for the exposed against 0.004 for the unexposed, and 19 of the 27 exposed in the top 27. Dependence on the outcome: in linear regression \(W = I\) and \(H = X(X'X)^{-1}X'\) is a function of the design alone. In a GLM \(w_i = \hat\mu_i(1 - \hat\mu_i)\) depends on the fitted mean, which depends on the coefficients, which depend on \(y\). So a patient’s leverage is also largest when their fitted risk is near 0.5 (where \(w\) peaks at 0.25) and shrinks toward the tails; the six top-leverage patients have fitted risks of 0.35 to 0.50. The same covariate row would carry less leverage in a cohort where it implied a risk of 0.02.

Q9. Flaw in the model, flaw in the data, or neither; what to report. Neither: it is a fact about the data that the model reports faithfully. The exposure coefficient is a function of the 6 exposed deaths and the 21 exposed survivors (the canonical-link identity of Q4 makes the fitted exposed deaths equal the observed six), so each exposed death carries about a sixth of the information about the target. Removing one takes the exposed death count from 6 to 5, and the OR from 1.54 to 1.22. Cook’s distance did not uncover an instability the printout hid: the standard error of 0.51 (a 95% interval of 0.57 to 4.16) already said that estimates of 1.22, or 2.0, are entirely compatible with these data, and p = 0.40 said the same thing in a third currency. Report the estimate with its interval, the number of exposed patients and exposed events, and the influence analysis as a description of where the information lives, not as a reason to delete anyone: the influential patients are the data the target is about, and deleting the six exposed deaths produces separation (the exposure coefficient runs to −14.9 with SE 487, because no exposed patient dies). The design-level cause is the one Lab 2 priced: 2.1% exposed. More exposed patients, not a different model, is the repair.

Understanding notes

  • The warnings are the lesson. The logbin-start chunk prints step size truncated due to divergence, step size truncated: out of bounds, and (for the 25-iteration call) glm.fit: algorithm did not converge. In the student render they appear as warning blocks under the chunk. Do not add #| warning: false to that chunk; the handout and the README tell students to expect them.
  • The boundary fit’s iteration count (29 with maxit = 500, 34 from the alternative start) is a property of R’s glm.fit step-halving on this platform (R 4.6.1). A future R version could change it; the estimate (RR 1.285, max fitted 0.9974) should not move. If the fit ever fails outright from the lecture’s start vector, the handout’s text needs rewording, not the start vector.
  • Poisson robust SE type. The handout uses HC3 (Lab 2’s convention and the one Week 7 recommends for small samples); Zou’s original modified Poisson uses HC0 (0.385 here). Either makes the point.
  • The rate model’s robust SE (0.532) is larger than its model-based one (0.426), the opposite direction from the no-offset Poisson. The reason is Lab 4’s Step 2 point: death ends follow-up, so expected deaths over a patient’s own follow-up are not Poisson-distributed around \(\text{rate} \times T\) in the way the model asserts. Key-only; the handout does not raise it.
  • The 1,271-row sample. Lab 3’s model dropped the three patients with no baseline heart rate (all unexposed survivors, ages 30, 43, 68), and this lab forms that sample explicitly so the four fits share rows. Consequence: Step 3’s numbers differ from the Lab 4 key’s (fit on 1,274) at the second decimal (deviance 820.25 vs 820.60, 815.67 vs 815.95, 759.76 vs 760.23; AICs likewise), while the by-age table rows quoted, the 73 ages, the 30 empty strata, the knots, and the 19 strata with an exposed patient are unchanged. The age-only exposure ORs of Lab 4’s Step 4 (2.163 / 2.137 / 2.375) become 2.160 / 2.135 / 2.375 on 1,271; the handout quotes the full-model versions (1.535 / 1.558 / 1.724) instead, so that the first row is Step 1’s logistic fit.
  • Canonical sums. colSums(preds) gives 136.00 for the logistic, Poisson, and LPM fits and 135.52 for the log-binomial. The Poisson sum is exact even though one fitted value exceeds 1. The rate model’s fitted values also sum to 136 (same canonical pair).
  • The deviance-residual plot’s two bands. For a 0/1 outcome, \(r_{D,i} = \sqrt{-2\log\hat\mu_i}\) for a death and \(-\sqrt{-2\log(1-\hat\mu_i)}\) for a survivor: two deterministic curves of the fitted risk. Students sometimes read the upper band’s slope as a trend; it is the formula.
  • The binned residual plot (Gelman & Hill 2007, §5.6) uses \(\pm 2\sqrt{\bar p(1-\bar p)/n}\) bands; with ten bins, about half a bin is expected outside under perfect calibration, so the second decile (0.071 observed vs 0.035 fitted, 2.2 SE) is not evidence of misfit by itself. The age-band version (key-only) shows the 30-to-39 band with zero deaths against 0.028 expected (1.4 SE).
  • hatvalues() uses fit$weights, the last iteration’s working weights, which differ from \(\hat\mu(1-\hat\mu)\) at the converged coefficients by ~1e−8. The handout’s by-hand version uses the same weights, hence the 1e−15 match. The rowSums((A %*% B) * A) idiom computes the diagonal of \(ABA'\) without forming the 1,271 × 1,271 matrix.
  • Cook’s distance for a GLM in R is \(r_{P,i}^2 / (p\,\hat\phi) \cdot h_i/(1-h_i)^2\) with \(\hat\phi = 1\) for the binomial, Pregibon’s (1981) one-step approximation to the change in \(\hat\beta\) measured in \((X'WX)\) units. The by-hand formula in the handout matches to 7e−18.
  • The 4/n rule flags 121 patients (every death with a low fitted risk has a Pearson residual of 3 or more), so it is useless for a rare binary outcome; the handout ranks instead of thresholding.
  • Deleting the top-1 patient moves the exposure coefficient by 0.23 and every other coefficient by at most 0.06: the influence is specific to the target, which is the point of Q9. Deleting the six exposed deaths gives separation (coefficient −14.9, SE 487, 20+ iterations); key-only.
  • The LPM residual-variance line (0.041 / 0.066 / 0.160 by tertile of fitted risk) is the forward pointer to Week 7; \(\mu(1-\mu)\) at the tertile means is 0.025 / 0.096 / 0.153, close for the upper two tertiles and above the model’s number in the lowest one (where the LPM’s fitted risks are partly negative). The handout only claims the variance grows.
  • Predicting a factor(age) fit works only at observed ages (the spectrum figure uses sat_grd); a continuous grid errors with “factor has new levels”.
  • Tibble printing. dat[top, ...] prints as a tibble; the handout avoids summarise() tables that round to three significant digits except in bins, which is wrapped in as.data.frame(round(...)).
  • Time. After the instructor’s 2026-10-13 trim (Step 0, about 4 minutes, and the original Q7) the step budget sums to about 60 minutes (3 + 16 + 10 + 15 + 16) against the 80-minute slot, with nine questions (Lab 3’s trim went 10 to 7; this one 10 to 9).