Lab 6 — Computing Variance

EPID 785R · Regression and Study Design

Author

YOUR NAME HERE

Published

October 25, 2026

Introduction

In the last lecture, we were introduced to variance estimation. The conceptual definition of a standard error is that it is an estimate of the standard deviation of the estimates we would get if we drew sample after sample from the same super-population and fit the same model to each.

The lecture notes closed with a table representing model-based, robust, and cluster-robust standard errors, all from the perspective of the “sandwich” variance estimator:

\[(X^{\intercal}X)^{-1} \, X^{\intercal}\Sigma X \, (X^{\intercal}X)^{-1}\]

What changes with each is that a different rule govern what fills \(\Sigma\).

We also the bootstrap, which replaces the analytic expression for the variance with resampling, treating the original data as a pseudo-source population.

In this lab we further explore the idea variance estimation from the sandwich (model-based and robust) and bootstrap perspectives. Unlike Labs 1 to 5, we won’t actually use the eICU data. Instead, we’ll simulate data.

This will allow us to generate actual data for the super-population, sample from it (or, at least, from the data generating mechanism that governs it), and explore how different variance estimators capture the true variability inherent in the sampling process.

We’ll see that the correct variance eestimators are the ones whose average over all the simulation samples matches the standard deviation of the estimates, and whose 95% intervals contain the truth 95% of the time. We’ll change two parameters governing the data generating mechanisms to get a sense of what each is doing: heteroscedasticity and sample size.

In this lab, we will:

draw one sample and compute the model-based and robust (HC0, HC3) standard errors of the exposure coefficient by hand, from the design matrix and the residuals, then check each against vcov() and sandwich::vcovHC(), see how the four HC variants differ, and reproduce what R, Stata, and SAS each print by default for the same model;
draw two thousand samples from worlds with and without heteroscedasticity, and compare each standard error’s average, and the coverage of its 95% interval, to the Monte Carlo truth;
compute the bootstrap standard error by hand and with the boot package, and put the bootstrap inside the simulation to check its coverage too.

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

Background for this lab: the Week 7 notes (Variance Estimation) and Mansournia et al. (2021) on the reading list. The by-hand sandwich in Step 1 follows the notes’ section 2 and 3; the bootstrap of Step 3 is section 4. Every simulation in this lab is the notes’ opening experiment (twenty thousand samples, one standard error) run again with a different standard error, or a different world.

1 Before you start (~3 min)

Open lab6.Rproj in RStudio.
There is no data to download. Every dataset in this lab is drawn by the function below.
Packages: the tidyverse, sandwich and lmtest (Labs 2 and 5), and boot (ships with R).

The function make_data() is our super-population. Read the function before running it to get a sense of how the data are generated. A confounder \(z\) is standard normal; the exposure \(x\) is more common at high \(z\) (about 21% exposed overall); the outcome \(y\) has mean \(1 + 2x + 1.5z\), so the truth we will try to recover is 2 is encoded in this outcome model as a mean difference of 2.

The error around that mean has a standard deviation that depends on exposure (heteroscedastic): 2 among the exposed, 1 among the unexposed, by default. The exposed are the rarer group and the more variable one, the mirror image of the lecture’s cluster trial, where the untreated practices were the more variable arm.

pacman::p_load(tidyverse, sandwich, lmtest, boot)

expit <- function(x) 1 / (1 + exp(-x))

# One sample from the super-population. The truth is written in the code:
# the exposure coefficient is 2, and the error SD depends on exposure.
make_data <- function(n, sd_exposed = 2, sd_unexposed = 1) {
  z     <- rnorm(n)                                    # a confounder
  x     <- rbinom(n, 1, expit(-1.5 + 0.8 * z))         # exposure: ~21%, more common at high z
  sigma <- ifelse(x == 1, sd_exposed, sd_unexposed)    # the error SD, by exposure group
  y     <- 1 + 2 * x + 1.5 * z + rnorm(n, 0, sigma)    # the outcome; beta_x = 2 is the truth
  data.frame(y, x, z)
}

2 Step 1 — One sample, the sandwich by hand (~17 min)

Draw one sample of 300 and fit the linear model we will use all morning (provided). We know the truth, so we can also check what the model cannot see: the residual standard deviation in each exposure group.

set.seed(785)
dat <- make_data(n = 300)
sum(dat$x)                                             # the exposed

fit <- lm(y ~ x + z, data = dat)
round(summary(fit)$coefficients, 4)
round(tapply(resid(fit), dat$x, sd), 3)                # residual SD, unexposed then exposed (truth 1 and 2)

Checkpoints: 63 exposed. The exposure coefficient is 2.313 with a model-based standard error of 0.193. The residual standard deviations are 1.07 among the unexposed and 1.97 among the exposed: again, heteroscedasticity (that we built in).

2.1 The model-based variance

Let’s start with the model-based variance, using the lecture’s approach, with three ingredients: the design matrix \(X\), the residuals \(\hat u = y - X\hat\beta\), and \(\Sigma\).

For the model-based estimator, \(\Sigma\) is \(s^2\) times the identity: one number on the diagonal for everyone. Complete TODO (1a): that number is the sum of the squared residuals divided by the residual degrees of freedom, \(N - p\).

X <- model.matrix(fit)            # n x 3: intercept, x, z
y <- dat$y
n <- nrow(X); p <- ncol(X)

beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
max(abs(beta_hat - coef(fit)))                          # zero: Week 4's formula

u     <- as.vector(y - X %*% beta_hat)                  # the residuals
bread <- solve(t(X) %*% X)                              # (X'X)^-1, the bread

# TODO (1a): the residual variance, one number for everyone
s2 <- sum(???) / (??? - ???)
s2
Sigma_model <- diag(s2, n)                              # s2 on the diagonal, zeros elsewhere

V_model <- bread %*% (t(X) %*% Sigma_model %*% X) %*% bread
max(abs(V_model - vcov(fit)))                           # zero
sqrt(V_model["x", "x"])                                 # the SE the printout showed

Checkpoints: \(s^2\) = 1.716, and the by-hand matrix equals vcov(fit) to machine precision, so the by-hand standard error of the exposure coefficient is the printout’s 0.1927. Note what \(s^2\) is: a pooled residual variance, about 1.7, sitting between the unexposed group’s 1.14 (= 1.07²) and the exposed group’s 3.81 (= 1.97²). It is too large for the one and too small for the other.

2.2 The robust variance

Keep the bread, change the meat. The HC0 estimator puts each person’s own squared residual on the diagonal of \(\Sigma\); HC3 divides each by \((1 - h_i)^2\), where \(h_i\) is the person’s leverage, the diagonal of the hat matrix \(X(X^{\intercal}X)^{-1}X^{\intercal}\) (Lab 5 computed the GLM version of this). Complete TODO (1b): the two diagonals.

h <- diag(X %*% bread %*% t(X))                         # leverage
c(max(abs(h - hatvalues(fit))), sum(h))                 # zero, and p = 3

# TODO (1b): the diagonal of Sigma under HC0 (each squared residual) and HC3
Sigma_hc0 <- diag(???, n)
Sigma_hc3 <- diag(??? / (1 - ???)^2, n)

V_hc0 <- bread %*% (t(X) %*% Sigma_hc0 %*% X) %*% bread
V_hc3 <- bread %*% (t(X) %*% Sigma_hc3 %*% X) %*% bread

max(abs(V_hc0 - vcovHC(fit, type = "HC0")))             # zero
max(abs(V_hc3 - vcovHC(fit, type = "HC3")))             # zero

se_tab <- c(model = sqrt(V_model["x", "x"]), HC0 = sqrt(V_hc0["x", "x"]), HC3 = sqrt(V_hc3["x", "x"]))
round(se_tab, 4)

round(range(u^2), 3)                                    # the HC0 diagonal, smallest to largest
round(tapply(u^2, dat$x, mean), 3)                      # its mean, unexposed then exposed

Checkpoints: both by-hand matrices match vcovHC() to machine precision. The standard error of the exposure coefficient is 0.193 model-based, 0.261 HC0, and 0.265 HC3. The HC0 diagonal runs from 0 to 25.2 (the model-based one was 1.716 for all 300 rows), and its average is 1.14 among the unexposed and 3.81 among the exposed: the robust meat has, person by person, recovered the two variances the model-based meat pooled away.

2.3 The HC family, and what the software returns

All five variants on this sample (provided): "const" is the model-based estimator, and the rest differ only in how each squared residual is rescaled before it enters \(\Sigma\) (per the notes’ table at the end: HC1 by \(N/(N-p)\), HC2 by \(1/(1-h_i)\), HC3 by \(1/(1-h_i)^2\)). The coeftest() line is a function we can use to print a coefficient table with a robust standard error (as we would if we looked at the coefficients from the model fit directly).

hc <- sapply(c("const", "HC0", "HC1", "HC2", "HC3"),
             function(type) sqrt(vcovHC(fit, type = type)["x", "x"]))
round(hc, 4)

round(range(h), 4); round(mean(h), 4); p / n              # leverage: range, mean, and p/n
round(tapply(h, dat$x, mean), 4)                        # mean leverage, unexposed then exposed

round(coeftest(fit, vcov = vcovHC(fit, type = "HC3"))["x", ], 4)

Checkpoints: 0.1927 / 0.2606 / 0.2619 / 0.2630 / 0.2654. The leverages run from 0.004 to 0.048 and average 0.010, which is \(p/n = 3/300\); the exposed carry two and a half times the leverage of the unexposed (0.019 against 0.008). With leverages this small, HC0 to HC3 are within two percent of one another, and all are a third larger than the model-based value.

Robust defaults differ across software

In R, sandwich::vcovHC() returns HC3 unless told otherwise, and sandwich::sandwich() returns HC0. Stata’s , vce(robust) is HC1. In SAS, the REPEATED statement of PROC GENMOD returns an HC0-type empirical variance, and PROC REG’s HCCMETHOD= option selects 0 through 3. The four agree as \(N\) grows (leverages shrink and \(N/(N-p)\) goes to one), and differ in small samples, where HC0 is known to be too small and HC3 is the recommended choice (Long & Ervin 2000). When two analysts report different “robust” standard errors from the same data, this is usually why.

2.4 The same model in Stata and SAS

Each default is a known formula, so we can print, from R, the standard error that each package would report for this model (provided). The chunk also writes the sample to lab6_sample.csv; anyone with Stata or SAS can run the code shown after the chunk and check the printout against the table.

write_csv(dat, "lab6_sample.csv")                      # the sample, for Stata or SAS

se_of <- function(type) sqrt(vcovHC(fit, type = type)["x", "x"])
tribble(
  ~software, ~command,                                      ~variant,
  "R",       "vcovHC(fit)                     (the default)", "HC3",
  "R",       "sandwich(fit)",                                 "HC0",
  "Stata",   "regress y x z, vce(robust)",                    "HC1",
  "Stata",   "regress y x z, vce(hc3)",                       "HC3",
  "SAS",     "proc reg; model y = x z / hcc hccmethod=0;",    "HC0",
  "SAS",     "proc reg; model y = x z / hcc hccmethod=3;",    "HC3",
  "SAS",     "proc genmod; ...; repeated subject=id;",        "HC0") |>
  mutate(SE = round(sapply(variant, se_of), 4)) |>
  as.data.frame()

The Stata and SAS code behind the table, for the record (not run here):

import delimited "lab6_sample.csv", clear
regress y x z                            // model-based:   0.1927
regress y x z, vce(robust)               // HC1:           0.2619
regress y x z, vce(hc3)                  // HC3:           0.2654

proc import datafile="lab6_sample.csv" out=lab6 dbms=csv replace; run;
proc reg data=lab6;
  model y = x z / hcc hccmethod=3;       * HC3: 0.2654 (hccmethod=0 gives HC0, 0.2606);
run;
data lab6; set lab6; id = _n_; run;      * one "subject" per row;
proc genmod data=lab6;
  class id;
  model y = x z / dist=normal;
  repeated subject=id;                   * empirical (sandwich) SE = HC0: 0.2606;
run;

Questions — write your answers below each prompt

Q1. Our by-hand model-based matrix matched vcov() and our by-hand HC3 matched vcovHC(). Looking at the two computations, name the one object that changed between them and the two that did not. Then: the model-based diagonal was 1.716 for everyone, while the residual variances were 1.14 among the unexposed and 3.81 among the exposed. For which group is the model-based \(\Sigma\) too small, and why would that matter more for the exposure coefficient than for the intercept? (Hint: which rows of \(X\) carry information about the exposure coefficient?)

YOUR ANSWER:

Q2. For one coefficient, we obtained five standard errors: 0.1927, 0.2606, 0.2619, 0.2630, 0.2654. Explain why HC0 < HC1 and HC0 < HC2 < HC3 is guaranteed in every sample, not just this one (hint: \(0 < h_i < 1\) and \(N/(N-p) > 1\)). A colleague who ran the same model in Stata with vce(robust) reports 0.262, another who used PROC GENMOD in SAS reports 0.261, and we report 0.265. Is any of us wrong? What about the sample would make the disagreement disappear, and what should each of us have written next to the word “robust”?

YOUR ANSWER:

3 Step 2 — Many samples: which standard error is right? (~18 min)

One sample told us the model-based and robust standard errors disagree. It cannot tell us which is right, because “right” means “matches the standard deviation of the estimates across repeated samples”, and we only have one estimate. So let’s draw two thousand and see what happens.

The engine (provided, with one TODO) has two parts. simulate() draws \(R\) samples of size \(n\) from make_data(), fits the model to each, and stores the estimate and four standard errors. summarise_sim() then reports, for each standard error, two things: its average across samples, to compare with the standard deviation of the estimates (the Monte Carlo truth); and its coverage, the share of samples in which the 95% Wald interval, estimate \(\pm 1.96 \times\) SE, contained the truth. Complete TODO (2a): the coverage indicator.

summarise_sim <- function(res, truth = 2) {
  se_cols <- grep("^se_", names(res), value = TRUE)
  out <- c(mean_estimate = mean(res$estimate), sd_estimate = sd(res$estimate))
  for (s in se_cols) {
    # TODO (2a): TRUE when the 95% Wald interval, estimate +/- 1.96 * se, contains the truth
    covered <- abs(??? - ???) <= 1.96 * res[[s]]
    out[paste0("mean_", s)]               <- mean(res[[s]])
    out[paste0("cov_", sub("se_", "", s))] <- mean(covered)
  }
  round(out, 3)
}

simulate <- function(R, n, sd_exposed, sd_unexposed, seed = 785) {
  set.seed(seed)
  res <- vector("list", R)
  for (r in 1:R) {
    d   <- make_data(n, sd_exposed, sd_unexposed)
    fit <- lm(y ~ x + z, data = d)
    res[[r]] <- c(estimate = coef(fit)[["x"]],
                  se_model = sqrt(vcov(fit)["x", "x"]),
                  se_hc0   = sqrt(vcovHC(fit, type = "HC0")["x", "x"]),
                  se_hc1   = sqrt(vcovHC(fit, type = "HC1")["x", "x"]),
                  se_hc3   = sqrt(vcovHC(fit, type = "HC3")["x", "x"]))
  }
  as.data.frame(do.call(rbind, res))
}

Four worlds (provided): no heteroscedasticity; the default world, in which the rarer exposed group is the more variable one; its mirror image, in which the common unexposed group is the more variable one; and the default world again at \(n = 80\). Each takes a few seconds to run.

sim_homo    <- simulate(R = 2000, n = 300, sd_exposed = 1, sd_unexposed = 1)
sim_hetero  <- simulate(R = 2000, n = 300, sd_exposed = 2, sd_unexposed = 1)
sim_reverse <- simulate(R = 2000, n = 300, sd_exposed = 1, sd_unexposed = 2)
sim_small   <- simulate(R = 2000, n = 80,  sd_exposed = 2, sd_unexposed = 1)

head(sim_hetero, 3)                                    # one row per sample; row 1 is Step 1's sample

scen <- rbind(homoscedastic           = summarise_sim(sim_homo),
              exposed_more_variable   = summarise_sim(sim_hetero),
              unexposed_more_variable = summarise_sim(sim_reverse),
              exposed_more_var_n80    = summarise_sim(sim_small))
t(scen)

Checkpoints: the first row of sim_hetero is Step 1’s sample (same seed): 2.313, 0.193, 0.261, 0.262, 0.265. Across the 2,000 samples, every mean estimate is within 0.01 of 2 (no standard error can fix bias, and here there is none to fix). Then, reading down the columns:

Homoscedastic: the estimates have SD 0.152; the model-based SE averages 0.149, HC0 0.148, HC3 0.150; coverage 0.95 for all of them. When the model’s variance assumption is true, every route arrives at the same place.
Exposed more variable: SD 0.272; the model-based SE averages 0.190 (70% of the truth) and covers only 83%; HC0 averages 0.264 and HC3 0.269, covering 94% and 95%.
Unexposed more variable: SD 0.203; the model-based SE now averages 0.275, a third too large, and covers 99%; HC3 averages 0.201 and covers 95%.
\(n = 80\): SD 0.540; model-based 0.370 (coverage 0.82); HC0 0.497 (0.92), HC1 0.507 (0.93), HC3 0.537 (0.94).

The figure (provided) is the lecture’s opening figure drawn in the world where the model is wrong: the histogram of the 2,000 estimates, with a normal curve whose spread is the SD of the estimates and two more whose spreads are the average model-based and HC3 standard errors.

ggplot(sim_hetero, aes(x = estimate)) +
  geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "grey85", color = "grey60") +
  stat_function(fun = dnorm, aes(color = "SD of the 2,000 estimates"),
                args = list(mean = mean(sim_hetero$estimate), sd = sd(sim_hetero$estimate)),
                linewidth = 1) +
  stat_function(fun = dnorm, aes(color = "Mean model-based SE"),
                args = list(mean = mean(sim_hetero$estimate), sd = mean(sim_hetero$se_model)),
                linewidth = 1, linetype = "dashed") +
  stat_function(fun = dnorm, aes(color = "Mean HC3 SE"),
                args = list(mean = mean(sim_hetero$estimate), sd = mean(sim_hetero$se_hc3)),
                linewidth = 1, linetype = "dotted") +
  geom_vline(xintercept = 2, color = "grey30") +
  scale_color_manual(values = c("SD of the 2,000 estimates" = "#0072B2",
                                "Mean model-based SE" = "#D55E00",
                                "Mean HC3 SE" = "#009E73")) +
  labs(x = "Estimated exposure coefficient (truth = 2)", y = "Density", color = NULL) +
  theme_classic() + theme(legend.position = "bottom")

In the lecture’s figure the two curves were indistinguishable. Here the dashed curve (model-based) is visibly too narrow for the histogram it claims to describe, and the dotted one (HC3) sits on top of the solid one.

Questions

Q3. In the heteroscedastic world the model-based SE averaged 0.190 against a true SD of 0.272 and covered 83%; in the mirror-image world it averaged 0.275 against 0.203 and covered 99%. Using what fills \(\Sigma\) in the model-based estimator (one pooled \(s^2\)) and who is rare in these data, explain why the model-based SE is too small in the first world and too large in the second. Why did the robust SE land near the truth in both worlds without being told which one it was in?

YOUR ANSWER:

Q4. The verdict in this step needed 2,000 samples and the true value of the coefficient. In practice we have one sample and no truth. In Step 1 the model-based and HC3 standard errors were 0.193 and 0.265. The lecture’s closing paragraph said a gap like that is “worth reading as a diagnostic of the model, and not only as a correction to its uncertainty”. What does the gap tell us about our model here? And what would it not tell us if the mean model, \(1 + 2x + 1.5z\), were the part we had wrong (hint: robust to what?)

YOUR ANSWER:

Q5. At \(n = 300\) the four HC variants covered within a point of one another; at \(n = 80\) HC0 covered 92% and HC3 94%, with HC1 (Stata’s default) between them. Which piece of the HC3 formula does the work in small samples, and what does it do to the squared residual of a person with high leverage? Why does the gap close as \(n\) grows? (Hint: the mean leverage is \(p/n\).)

YOUR ANSWER:

4 Step 3 — The bootstrap (~15 min)

The bootstrap is an alternative route to variance estimation, one that does not require analytic expressions for the variance. Instead of filling \(\Sigma\) by some mathematical rule, it re-runs the sampling, with the observed sample standing in for the super-population: draw \(N\) rows with replacement, refit, keep the coefficient, repeat \(B\) times, and take the standard deviation of the \(B\) replicates. We do it by hand on Step 1’s sample first. Complete TODO (3a): the resampled row numbers.

set.seed(2026)
B <- 1000
b_star <- numeric(B)
for (b in 1:B) {
  # TODO (3a): n row numbers, drawn with replacement
  idx <- sample(???, replace = ???)
  b_star[b] <- coef(lm(y ~ x + z, data = dat[idx, ]))[["x"]]
}

round(sd(b_star), 4)                                    # the bootstrap standard error
round(se_tab, 4)                                        # model, HC0, HC3 from Step 1
round(mean(b_star) - coef(fit)[["x"]], 4)               # the bootstrap estimate of bias

round(coef(fit)[["x"]] + c(-1, 1) * 1.96 * sd(b_star), 4)   # normal-interval bootstrap CI
round(quantile(b_star, c(0.025, 0.975)), 4)                 # percentile bootstrap CI
round(coef(fit)[["x"]] + c(-1, 1) * 1.96 * se_tab["model"], 4)   # Wald, model-based
round(coef(fit)[["x"]] + c(-1, 1) * 1.96 * se_tab["HC3"], 4)     # Wald, HC3

Checkpoints: the bootstrap standard error is 0.275, against 0.193 model-based and 0.265 HC3. The bias estimate is 0.015. The normal-interval bootstrap CI is (1.78, 2.85) and the percentile CI (1.80, 2.85); the Wald intervals are (1.94, 2.69) with the model-based SE and (1.79, 2.83) with HC3. Notice, nowhere in the bootstrap loop did we mention heteroscedasticity, but the bootstrap landed on the robust variance estimate anyway.

The boot package does the same with a function of the data and an index vector (provided; the lecture’s boot_func). Its normal interval is centered on the bias-corrected estimate, which is why it differs slightly from the one we computed by hand; a second difference is that this run drew its own 1,000 resamples.

boot_stat <- function(data, index) {
  coef(lm(y ~ x + z, data = data[index, ]))[["x"]]
}

set.seed(2026)
boot_res <- boot(data = dat, statistic = boot_stat, R = 1000)
boot_res
boot.ci(boot_res, type = c("norm", "perc", "bca"))

Checkpoints: standard error 0.267, bias 0.013; normal (1.776, 2.824), percentile (1.821, 2.843), BC\(_a\) (1.812, 2.837). The two bootstrap runs (ours and the package’s) differ by 0.008 in the standard error: that is Monte Carlo error, with 1,000 replicates, and it shrinks as \(B\) grows.

4.1 The bootstrap inside the simulation

Last, the check we ran on the sandwich estimators in Step 2, run on the bootstrap (provided): 500 samples from the heteroscedastic world, 200 resamples each, so 100,000 fits. To keep that under a few seconds, the inner loop uses the Week 4 formula, \((X^{\intercal}X)^{-1}X^{\intercal}y\), directly instead of lm().

fast_ols <- function(X, y) solve(crossprod(X), crossprod(X, y))   # (X'X)^-1 X'y, Week 4

simulate_boot <- function(R, n, B, sd_exposed = 2, sd_unexposed = 1, seed = 785) {
  set.seed(seed)
  res <- vector("list", R)
  for (r in 1:R) {
    d <- make_data(n, sd_exposed, sd_unexposed)
    X <- cbind(1, d$x, d$z)
    b_hat  <- fast_ols(X, d$y)[2]
    b_star <- numeric(B)
    for (k in 1:B) {
      idx <- sample(n, replace = TRUE)
      b_star[k] <- fast_ols(X[idx, ], d$y[idx])[2]
    }
    fit <- lm(y ~ x + z, data = d)
    res[[r]] <- c(estimate   = b_hat,
                  se_model   = sqrt(vcov(fit)["x", "x"]),
                  se_hc3     = sqrt(vcovHC(fit, type = "HC3")["x", "x"]),
                  se_boot    = sd(b_star),
                  perc_lower = quantile(b_star, 0.025, names = FALSE),
                  perc_upper = quantile(b_star, 0.975, names = FALSE))
  }
  as.data.frame(do.call(rbind, res))
}

sim_boot <- simulate_boot(R = 500, n = 300, B = 200)
summarise_sim(sim_boot)
mean(sim_boot$perc_lower <= 2 & 2 <= sim_boot$perc_upper)   # percentile-interval coverage

Checkpoints: SD of the estimates 0.276; the model-based SE averages 0.190 (coverage 0.83), HC3 0.271 (0.96), and the bootstrap 0.267 (0.95); the percentile interval covers 0.93. (With 500 samples the Monte Carlo standard error of a coverage estimate is about 0.01, so the last three are not distinguishable from one another.)

Questions

Q6. The bootstrap SE on Step 1’s sample (0.275) matched the robust SE (0.265), not the model-based one (0.193), and in the simulation it covered 95% in a world where the model-based SE covered 83%. We never told it the exposed were more variable. Why does it not need to be told? Then: there is one thing the loop did have to be told, in the line we completed. What is it, and what would go wrong if the rows of dat were not independent draws?

YOUR ANSWER:

Q7. In this lab the normal-interval and percentile bootstrap intervals nearly coincide. The lecture said the percentile interval is transformation- and range-respecting and the normal interval is not. Name one quantity we have estimated in Labs 2 to 5 for which that distinction would matter, and say what we would do about it with the normal interval. Then a caution: the Lab 1 cohort has 27 exposed patients, 6 of whom died. If we bootstrapped the exposure log odds ratio there, what would some of the resamples contain, and what would that do to the replicates? Which of the three routes to a standard error is most exposed to that problem, and is any of them safe?

YOUR ANSWER:

5 Bonus — Lab 5’s question, answered by simulation (~5 min)

Lab 5’s Q3 asked why the Poisson model’s model-based standard error for a binary outcome was 60% too large in the lecture’s NHEFS example and only 5% too large in the ICU cohort, where the risk of death was 0.11. Here is that question as a simulation (provided): a binary outcome with a log-linear risk (true risk ratio 1.3), fit with the Poisson family and the log link (Zou’s route to the risk ratio, Lab 2 and Lab 5), in a world where the risk is about 0.55 and in one where it is about 0.11. The robust SE is HC0, the one the lecture derived by hand for a GLM.

make_binary <- function(n, risk_unexposed, rr = 1.3) {
  z  <- rnorm(n)
  x  <- rbinom(n, 1, expit(-0.85 + 0.5 * z))               # about 30% exposed
  mu <- pmin(risk_unexposed * rr^x * exp(0.2 * z), 0.98)   # a log-linear risk; RR = 1.3
  y  <- rbinom(n, 1, mu)
  data.frame(y, x, z)
}

simulate_poisson <- function(R, n, risk_unexposed, seed = 785) {
  set.seed(seed)
  res <- vector("list", R)
  for (r in 1:R) {
    d   <- make_binary(n, risk_unexposed)
    fit <- glm(y ~ x + z, family = poisson(link = "log"), data = d)
    res[[r]] <- c(estimate  = coef(fit)[["x"]],
                  se_model  = sqrt(vcov(fit)["x", "x"]),
                  se_robust = sqrt(vcovHC(fit, type = "HC0")["x", "x"]),
                  mean_risk = mean(d$y))
  }
  as.data.frame(do.call(rbind, res))
}

pois_common <- simulate_poisson(R = 1000, n = 500, risk_unexposed = 0.5)
pois_rare   <- simulate_poisson(R = 1000, n = 500, risk_unexposed = 0.1)

t(rbind(risk_about_0.55 = summarise_sim(pois_common, truth = log(1.3)),
        risk_about_0.11 = summarise_sim(pois_rare,   truth = log(1.3))))
round(c(mean(pois_common$mean_risk), mean(pois_rare$mean_risk)), 3)

Checkpoints: both mean estimates are within 0.01 of \(\log(1.3) = 0.262\). At a risk of 0.56 the estimates have SD 0.079, the model-based SE averages 0.127 (60% too large; coverage 1.000) and the robust SE 0.077 (0.94). At a risk of 0.11 the SD is 0.270, the model-based SE averages 0.285 (5% too large; coverage 0.97) and the robust SE 0.266 (0.95). (One more software note, since this is the model most often run in Stata as poisson y x z, vce(robust): for its maximum-likelihood commands Stata multiplies HC0 by \(N/(N-1)\), not \(N/(N-p)\), so its “robust” standard error for a GLM is HC0 times 1.001 at \(N = 500\), and sandwich::vcovHC(fit, type = "HC1") does not reproduce it exactly; vcovHC(fit, type = "HC0") * N/(N-1) does.)

Questions

Q8. Write the answer to Lab 5’s Q3 now, with the simulation in hand: the model-based meat for a Poisson model holds \(\hat\mu_i\) on its diagonal, the true variance of a binary outcome is \(\mu_i(1-\mu_i)\). What is the ratio of the two at \(\mu = 0.56\) and at \(\mu = 0.11\), and how does the square root of that ratio compare to the model-based SE’s excess in each world? In which direction does the model-based SE err for a Poisson model of a binary outcome, always, and is that the same direction as the error in Step 2?

YOUR ANSWER:

References

Mansournia MA, Nazemipour M, Naimi AI, Collins GS, Campbell MJ. Reflection on Modern Methods: Demystifying Robust Standard Errors for Epidemiologists. Int J Epidemiol. 2021;50(1):346–351. (Reading list; the by-hand sandwich of Step 1.)
White H. A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity. Econometrica. 1980;48(4):817–838. (HC0.)
MacKinnon JG, White H. Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties. J Econometrics. 1985;29(3):305–325. (HC1 to HC3.)
Long JS, Ervin LH. Using Heteroscedasticity Consistent Standard Errors in the Linear Regression Model. Am Stat. 2000;54(3):217–224. (Why HC3 in small samples, Step 2.)
Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Chapman & Hall, 1993. (Step 3.)
Zeileis A. Econometric Computing with HC and HAC Covariance Matrix Estimators. J Stat Softw. 2004;11(10):1–17. (The sandwich package’s vcovHC().)
Zou G. A Modified Poisson Regression Approach to Prospective Studies with Binary Data. Am J Epidemiol. 2004;159(7):702–706. (The Bonus.)
Freedman DA. On the So-Called “Huber Sandwich Estimator” and “Robust Standard Errors”. Am Stat. 2006;60(4):299–302. (Q4.)