Lab 6 — Instructor answer key

Companion file: solution.R (in this folder), the lab chunk code with the four TODOs (1a, 1b, 2a, 3a) 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 (ratios of mean SE to Monte Carlo SD, the implied coverage of an SE that is 70% of the truth, and 1/√(1−μ) for the Bonus). About 20 seconds in all.

Contents: TODO code answers · expected checkpoints · model answers Q1–Q8 · understanding notes

TODO code answers

TODO	Answer
1a	`s2 <- sum(u^2) / (n - p)`, the residual sum of squares over the residual degrees of freedom (300 − 3). Gives 1.7162. Students who divide by `n` get 1.699 and a by-hand matrix that misses `vcov(fit)` by about 1%; the lecture’s `s` used `nrow(bmi_data) - 2` for the same reason. `summary(fit)$sigma^2` is the same number.
1b	`Sigma_hc0 <- diag(u^2, n)` and `Sigma_hc3 <- diag(u^2 / (1 - h)^2, n)`. Both by-hand matrices match `vcovHC()` to 1e−15. Students who write `diag(u^2)` (no `n`) get the same matrix (R builds a diagonal matrix from a vector of length n); students who write `diag(u, n)^2` get the same diagonal too but should be shown that squaring a matrix elementwise is not the matrix square.
2a	`covered <- abs(res$estimate - truth) <= 1.96 * res[[s]]`. Equivalent: `res$estimate - 1.96 * res[[s]] <= truth & truth <= res$estimate + 1.96 * res[[s]]`. The `truth` argument defaults to 2 and is set to `log(1.3)` in the Bonus. A student who hard-codes `2` gets the Bonus table wrong (coverage near zero), which is a useful mistake to see.
3a	`idx <- sample(nrow(dat), replace = TRUE)`, or `sample(1:300, 300, replace = TRUE)` or `sample(n, replace = TRUE)` (`n` is 300 from Step 1). With `replace = FALSE` every resample is the original sample in a different order and the SD of the replicates is 0.

Expected checkpoints

Verified 2026-10-23 on R 4.6.1. Every random chunk is seeded (785 for samples and simulations, 2026 for bootstraps), so these reproduce exactly under R’s default RNG.

setup:   prevalence in 100,000 draws 0.207; error SD by group 0.997 / 2.003
Step 1:  one sample n = 300: 63 exposed; coefficients 1.0002 / 2.3134 / 1.4464, model SEs 0.0851 / 0.1927 / 0.0782
         residual SD by group 1.070 (unexposed) / 1.966 (exposed)
         by hand: beta_hat = coef(fit) to 1e-14; s2 = 1.7162; V_model = vcov(fit) to 2e-16; SE(x) 0.1927
         leverage = hatvalues to 7e-16, sums to 3; range 0.0042-0.0475, mean 0.0100 = p/n; by group 0.0076 / 0.0192
         HC0 and HC3 by hand = vcovHC to 8e-16; SE(x) model 0.1927 / HC0 0.2606 / HC3 0.2654
         HC0 diagonal 0 to 25.166; its mean by group 1.139 / 3.806 (= the group residual variances, 1.07^2 and 1.97^2)
         HC family: const 0.1927 / HC0 0.2606 / HC1 0.2619 / HC2 0.2630 / HC3 0.2654; HC1 = HC0 * sqrt(300/297)
         coeftest HC3: estimate 2.3134, SE 0.2654, t 8.72, p 2e-16
         software defaults (by formula; writes lab6_sample.csv): R vcovHC 0.2654 (HC3) / sandwich() 0.2606 (HC0) /
           Stata regress vce(robust) 0.2619 (HC1) / Stata vce(hc3) 0.2654 / SAS REG hccmethod=0 0.2606, =3 0.2654 /
           SAS GENMOD repeated subject=id 0.2606 (HC0)
         key-only: sandwich(fit) = HC0 exactly; vcovCL with one row per cluster and cadjust = FALSE = HC0 (the GENMOD device);
           had this been a Stata ML command (glm/poisson), vce(robust) = HC0 * sqrt(n/(n-1)) = 0.2611
Step 2:  four scenarios, 12 s; row 1 of sim_hetero = Step 1's sample (2.3134, 0.1927, 0.2606, 0.2619, 0.2654)
         scenario                 mean_est  sd_est  se_model(cov)  se_hc0(cov)   se_hc1(cov)   se_hc3(cov)
         homoscedastic            2.001     0.152   0.149 (0.953)  0.148 (0.949) 0.149 (0.951) 0.150 (0.953)
         exposed_more_variable    2.002     0.272   0.190 (0.830)  0.264 (0.944) 0.265 (0.945) 0.269 (0.950)
         unexposed_more_variable  2.001     0.203   0.275 (0.990)  0.198 (0.946) 0.199 (0.946) 0.201 (0.949)
         exposed_more_var_n80     2.007     0.540   0.370 (0.824)  0.497 (0.919) 0.507 (0.925) 0.537 (0.936)
         key-only: se_model / sd_est = 0.98 / 0.70 / 1.36 / 0.69; se_hc3 / sd_est = 0.99 in all four
         key-only: an SE at 70% of the truth implies coverage P(|Z| < 1.96 * 0.70) = 0.83
         key-only: exposed count across the 2,000 n = 80 samples: 5 to 28 (never 0 or 1, so HC3 never divides by zero)
Step 3:  by-hand bootstrap, B = 1000, seed 2026: SD 0.2745 (model 0.1927 / HC0 0.2606 / HC3 0.2654); bias 0.0147
         normal (1.7753, 2.8515); percentile (1.7953, 2.8528); Wald model (1.9358, 2.6911); Wald HC3 (1.7932, 2.8336)
         boot(): SE 0.2673, bias 0.0132; normal (1.776, 2.824); percentile (1.821, 2.843); BCa (1.812, 2.837)
         boot sim 500 x 200, 2 s: sd 0.276; model 0.190 (0.828); HC3 0.271 (0.956); boot 0.267 (0.952); percentile 0.932
         key-only: MC SE of a coverage estimate: 0.010 at R = 500, 0.005 at R = 2000
Bonus:   two scenarios, 3 s; truth log(1.3) = 0.2624
         risk_about_0.55: mean_est 0.260, sd 0.079, model 0.127 (1.000), robust 0.077 (0.940); mean risk 0.558
         risk_about_0.11: mean_est 0.251, sd 0.270, model 0.285 (0.966), robust 0.266 (0.953); mean risk 0.113
         key-only: 1/sqrt(1 - mu) = 1.505 / 1.062 vs se_model / sd_est = 1.608 / 1.056

Model answers

Any phrasing that names the right mechanism is fine. Students should quote the numbers; the point of a simulation lab is that the numbers are the argument.

Q1. The one object that changed is $\Sigma$, the meat: $s^2 I$ for the model-based estimator, $\mathrm{diag}(\hat u_i^2 / (1-h_i)^2)$ for HC3. The two that did not change are the bread, $(X^{\intercal}X)^{-1}$, and the coefficient vector $\hat\beta$ (and therefore the residuals themselves): a robust standard error changes what we say about a coefficient’s precision, not the coefficient. The model-based diagonal, 1.716, is too small for the exposed (whose residual variance is 3.81) and too large for the unexposed (1.14). It matters more for the exposure coefficient because the exposure column of $X$ is nonzero only in the 63 exposed rows: the exposure coefficient is, essentially, the adjusted difference between the exposed mean and the unexposed mean, and the variance of a difference in means is $\sigma_1^2/n_1 + \sigma_0^2/n_0$, dominated by the smaller group’s term. The model-based estimator replaces $\sigma_1^2 = 3.8$ in that dominant term by 1.7. The intercept is anchored by the 237 unexposed rows, for which 1.7 is too large, so the model-based SE of the intercept errs in the other direction (in this sample: 0.0851 model-based against 0.0697 HC3; the confounder’s SE moves the same way, 0.0782 to 0.0739).

Q2. HC1 $= $ HC0 $\times N/(N-p)$ and $N/(N-p) > 1$, so HC1 > HC0 in every sample. HC2 divides each $\hat u_i^2$ by $(1 - h_i)$ and HC3 by $(1-h_i)^2$; because every leverage satisfies $0 < h_i < 1$ (the hat matrix is a projection; $h_i$ is its diagonal, between $1/N$ and 1), each HC2 term exceeds the HC0 term and each HC3 term exceeds the HC2 term, and a sandwich with a larger meat (in the positive-definite sense) gives a larger variance for every coefficient. So HC0 < HC2 < HC3 always. HC1 is not ordered against HC2 and HC3 in general: when every leverage equals $p/N$ exactly, as in the lecture’s balanced cluster trial, HC1 = HC2; here HC1 (0.2619) sits between HC0 and HC2. None of the three is wrong. Each reports a consistent estimator of the same quantity with a different finite-sample rescaling: the Stata colleague’s 0.262 is HC1 (regress, vce(robust) multiplies HC0 by $N/(N-k)$), the SAS colleague’s 0.261 is HC0 (PROC GENMOD’s empirical sandwich carries no correction), and our 0.265 is HC3 (R’s vcovHC() default). Each of us should have written the variant next to the word “robust” (a “robust SE” is not fully specified without it), and ideally the software and command, since the same word maps to different formulas in each. What would make the disagreement disappear is a larger sample: as $N$ grows the leverages shrink toward zero (their mean is $p/N$, already 0.01 here) and $N/(N-p)$ toward one, so all four converge to HC0. In this sample the spread is already under two percent; at $n = 80$ (Step 2) it is eight percent.

Q3. The model-based meat uses one pooled residual variance, $s^2 \approx n_0\sigma_0^2/N + n_1\sigma_1^2/N$, a size-weighted average that leans toward the larger group (here 0.79 × 1 + 0.21 × 4 ≈ 1.6 to 1.7). The sampling variance of the exposure coefficient, $\sigma_1^2/n_1 + \sigma_0^2/n_0$, leans toward the smaller group, because $1/n_1$ is the big weight. In the default world the smaller group is the more variable one: the true variance is dominated by $\sigma_1^2 = 4$, and the model-based estimator uses 1.7 in its place, so it is too small (0.190 against 0.272; coverage 0.83, which is exactly $P(|Z| < 1.96 \times 0.70)$). In the mirror image the smaller group has $\sigma_1^2 = 1$ while $s^2 \approx 0.79 \times 4 + 0.21 \times 1 ≈ 3.4$, so the model-based estimator is too large (0.275 against 0.203; coverage 0.99). The rule: the model-based SE follows the variance of the common group; the truth follows the variance of the rare group. The robust meat uses each person’s own $\hat u_i^2$, whose expectation is (close to) that person’s own $\sigma_i^2$; so $X^{\intercal}\Sigma X$ is approximately $\sum_i x_i x_i^{\intercal}\sigma_i^2$ whichever way the variances fall, and the estimator does not need to be told. (The lecture’s balanced-trial remark is the special case where the two weightings coincide: with $n_1 = n_0$ the pooled $s^2$ gives the slope’s variance exactly right even under heteroscedasticity. Our design is unbalanced, which is why it fails here.)

Q4. A model-based SE and a robust SE estimate the same thing under the model’s assumptions, so a 40% gap (0.193 against 0.265) is evidence that the model’s variance assumption, constant residual variance, is false in these data, and a pointer to where to look: in Step 1 the residual variance by exposure group (1.14 vs 3.81) found it at once, as Lab 5’s tertile check did for the linear probability model. In that sense the gap is a diagnostic: it told us the model is misspecified in its variance part before any simulation did. What it does not tell us is that the mean model is right. The robust SE is robust to the variance being wrong, not to the mean being wrong (the lecture’s closing caution; Freedman 2006, King & Roberts 2015). If $E(y \mid x, z)$ were not $1 + 2x + 1.5z$, the coefficient would converge to the best linear approximation (Week 4’s projection), and the robust SE would be the honest SE for that quantity, not for the one we had in mind. Worse, a wrong mean model itself produces residuals whose spread varies with the covariates, so a large model-based/robust gap can be a symptom of mean misspecification as well as of heteroscedasticity. Either way, the response to a large gap is to look at the model, not only to swap the standard error.

Q5. The $1/(1-h_i)^2$ factor. Least-squares residuals are smaller than the errors they estimate, and more so at high leverage: under constant variance $\mathrm{Var}(\hat u_i) = \sigma^2(1 - h_i)$, because a high-leverage observation pulls its own fitted value toward itself. HC0 therefore underestimates the meat, most for the points that matter most. HC2 undoes the shrinkage exactly (dividing by $1 - h_i$ restores $\sigma^2$ in expectation under homoscedasticity); HC3 divides by $(1-h_i)^2$, overcorrecting slightly, which approximates the jackknife and is what makes it the recommended small-sample choice. At $n = 80$ with $p = 3$ the mean leverage is 0.0375 and an exposed person’s is roughly twice that, so $(1-h_i)^{-2}$ is 1.08 to 1.17 for the rows that drive the exposure coefficient: HC3 averages 0.537 against HC0’s 0.497 and covers 0.94 against 0.92. At $n = 300$ the mean leverage is 0.01, the factors are within 2–4% of one, and the variants coincide (and $N/(N-p)$, HC1’s factor, is 1.01). All of them converge to HC0 as $n$ grows because the leverages go to zero at rate $p/n$.

Q6. The bootstrap resamples whole rows, $(y_i, x_i, z_i)$ together. The resampling distribution of $\hat\beta$ therefore reflects whatever the joint distribution of the rows actually contains, including the larger spread of $y$ among the exposed. It never computes a pooled $s^2$ and never assumes one; in large samples the variance of the pairs-bootstrap replicates is the sandwich with squared residuals in the meat, which is why 0.275 landed next to HC3’s 0.265 and HC0’s 0.261 rather than the model-based 0.193, and why its coverage in the simulation (0.95) matched HC3’s. The one thing we did tell it is the resampling unit: sample(nrow(dat), replace = TRUE) says that the row is what was drawn independently from the super-population, so rows are what get redrawn. That is the bootstrap’s only assumption, and it is an assumption about the design, not about the variance. If the rows were not independent draws (the lecture’s cluster trial, in which patients in the same practice share a deviation), resampling rows would treat two patients from the same practice as independent draws, the dependence that inflates the true variance would be destroyed in the bootstrap world, and the replicates would spread like independent data, reproducing the model-based and HC answer rather than the cluster-robust one (the lecture’s section 3.2). The repair is to resample the unit that was sampled independently in the design, here the practice, keeping each practice’s rows together. The bootstrap reproduces whatever dependence its resampling preserves and destroys whatever its resampling breaks.

Q7. A risk ratio or odds ratio (Lab 2’s case-cohort RR, Lab 5’s OR of 1.54 and RRs of 1.29 and 1.41). The sampling distribution of a ratio is right-skewed and bounded below by zero; a normal interval computed from the SD of the ratio replicates is symmetric about the point estimate and can reach below zero, and $\exp(\mathrm{SD}(\log \hat\theta^*))$ is not $\mathrm{SD}(\hat\theta^*)$. With the normal interval we would do what the lecture does: compute the SD of the log-ratio replicates, form the interval on the log scale, and exponentiate its endpoints. The percentile interval gives the same answer either way, because the 2.5th percentile of $\exp(\theta^*)$ is the exponential of the 2.5th percentile of $\theta^*$. The caution: with 27 exposed and 6 exposed deaths in 1,271 rows, a resample’s number of exposed deaths is roughly Poisson with mean 6, so about one resample in 400 contains no exposed deaths (and more contain one or two), and the resampled count of exposed people ranges from about 15 to 40. Replicates with no exposed deaths have a log odds ratio of $-\infty$ (R returns a coefficient near −17 with a huge SE: Lab 5’s separation); the others are heavy-tailed. The SD of the replicates, and so the normal-interval bootstrap, is driven by those few resamples and is unstable; the percentile interval is less affected but still rests on a handful of events in its tails. The bootstrap is the route most exposed to this, because it re-enacts the sparseness. The robust sandwich with 27 exposed is not safe either: its meat sums 27 informative score contributions, and HC0 is known to be too small in such samples. The model-based SE from the logistic model is, for once, the most stable of the three here, because the binomial variance function is the right one for a binary outcome and it borrows strength from the model rather than from the six deaths; but all three rely on large-sample arguments that six events do not support. None is “safe”; the honest tools at that point are exact or profile-likelihood intervals, or a penalized fit (Firth, Week 10), and a sentence in the report saying so. Lab 3’s p = 0.40 and SE of 0.51 were already telling us this.

Q8. The ratio of the Poisson meat to the true one is $\hat\mu_i / [\mu_i(1-\mu_i)] = 1/(1-\mu_i)$. At $\mu = 0.56$ that is 2.27 with square root 1.51; at $\mu = 0.11$ it is 1.12 with square root 1.06. The simulation’s excess is 0.127 / 0.079 = 1.61 and 0.285 / 0.270 = 1.06: the second matches the formula exactly, the first is a little larger because the exposed in the common-outcome world have risks up to the 0.98 cap, where $1/(1-\mu)$ is enormous and the person-weighted ratio exceeds the ratio at the mean risk. The direction is always the same: $\mu > \mu(1-\mu)$ for every $\mu$ in (0, 1), so the model-based Poisson SE for a binary outcome is always too large (conservative: coverage 1.000 and 0.966), and by an amount that vanishes as the outcome becomes rare. That is Lab 5’s answer: the ICU cohort’s risk was 0.11, so the Poisson SE was only 5% too large there, against 60% at the NHEFS risk of 0.5. This differs from Step 2, where the direction depended on which group was more variable and could go either way. The robust SE (HC0, the lecture’s by-hand GLM sandwich) is right at both risks (0.940, 0.953), which is Zou’s (2004) argument for pairing the Poisson model with it.

Understanding notes

Seeds and Monte Carlo error. Every random chunk calls set.seed() (785 for samples and simulations, 2026 for bootstraps), so the checkpoints reproduce exactly under R’s default RNG (any R ≥ 3.6). If a student’s numbers differ in the third decimal, check RNGkind(); if they differ in the first, a seed was skipped or a chunk run twice. The simulations are estimates: the Monte Carlo SE of a coverage near 0.95 is 0.005 at R = 2,000, 0.007 at R = 1,000, and 0.010 at R = 500, and the relative SE of an estimated SD is about $1/\sqrt{2R}$ (1.6% at R = 2,000). So “HC3 covers 0.950, HC0 0.944” in the heteroscedastic world is not a real difference; “0.936 vs 0.919” at n = 80 is borderline on its own but the two share the same samples, which makes the paired difference more reliable than the marginal SEs suggest.
The by-hand bootstrap and boot() disagree slightly (0.2745 vs 0.2673). Same data, same seed, different resamples, because boot() draws its indices differently. Each is within one Monte Carlo SE (about 0.006 at B = 1,000) of the other. A student who expects identical numbers has learned something about the bootstrap.
HC3 and leverage one. If a group has a single member, that row has $h_i = 1$ and HC3 divides by zero (NaN). Across the 2,000 samples of size 80 the exposed count ranges from 5 to 28, so this never happens in the handout; a student who tries n = 40 will meet it (as the prototype did). It is a real limitation worth mentioning if asked: HC3 needs every covariate pattern to have company.
The balanced-design subtlety. The lecture’s cluster trial showed heteroscedasticity aligned with a balanced binary covariate leaving the slope’s model-based SE correct (pooled $s^2 \times 2/n$ equals $\sigma_1^2/n + \sigma_0^2/n$ exactly when $n_1 = n_0$). Our 21% design is deliberately unbalanced so that the model-based SE fails. A student who edits make_data() to 50% exposed will see the model-based coverage climb back to about 0.93 (not 0.95, because the confounder z and the sample-to-sample imbalance keep it from being exact); this is a good extension question, not a bug.
replicate() and .... The engine uses explicit for loops. replicate(R, one_rep(n, ...)) inside a function silently misreads the dots (a known base-R gotcha) and was the cause of a prototype bug; do not “simplify” the loops into replicate().
Runtime. Step 2 is the slow part (four vcovHC() calls per fit, 8,000 fits, 12 s); Step 3 and the Bonus are 2 to 3 s each. If a classroom laptop struggles, drop R to 1,000 in Step 2; the lessons survive, the third decimals do not.
The cross-software table is reproduced by formula, not run. Neither Stata nor SAS was available when the lab was built, so the table uses vcovHC() options that implement each package’s documented default: Stata regress, vce(robust) = HC1 ($N/(N-k)$ times HC0) and vce(hc2), vce(hc3) by name; Stata’s ML commands (glm, poisson, logit) with vce(robust) = HC0 times $N/(N-1)$, which vcovHC(type = "HC1") does not reproduce (it uses $N/(N-p)$); SAS PROC REG HCCMETHOD=0..3 = HC0..HC3; SAS PROC GENMOD REPEATED SUBJECT=id with one subject per row = HC0 (verified with vcovCL(cluster = seq_len(n), type = "HC0", cadjust = FALSE)). The CSV the handout writes (lab6_sample.csv) lets anyone with either package confirm the printout; the Stata and SAS snippets in the handout are written against that file. If a student with Stata reports a different number, the first thing to check is whether they used regress or glm (the finite-sample factors differ). GENMOD estimates a scale parameter for dist=normal, but the empirical SE does not depend on it.
What this lab does not do. It no longer touches clustering: the cluster trial, the cluster-robust sandwich by hand and with vcovCL(), the cluster simulation, the cluster bootstrap, and the old Q6–Q7 were cut by the instructor on 2026-10-25 for length and are kept, with code, checkpoints, and model answers, in retired_step3_clusters.md in this folder (never ships). It does not fit a GEE (the lecture did; geepack is not loaded), does not show the parametric bootstrap, and does not compute BC$_a$ by hand. It also never shows a biased estimator: every mean_estimate is 2.00, so that the whole lab is about the second-order problem and students cannot confuse an SE repair with a bias repair. Q4 and the Introduction carry that message in words.