Same Coefficient, Different Standard Error: The Sandwich Estimator Worked by Hand

อ่านฉบับภาษาไทย (Thai version)

Abstract

Estimating a regression coefficient and estimating its uncertainty are two separate processes. A robust standard error is not a new model; it is a different estimate of the variance-covariance matrix under the same Poisson model. In a hypothetical study of treatment and hospital admissions, the coefficient stays at 0.40 while the standard error rises from 0.10 to 0.19. Two datasets of four patients, both with a mean of 5 admissions, show why. Poisson regression gives both the coefficient log 5 = 1.609 and the same conventional standard error, 0.224, because it assumes that the variance equals the mean. The robust (HC0) standard error reads the squared residuals: 0.071 for the tightly packed counts and 0.400 for the spread-out ones. The sandwich formula, bread and meat, extends this to many predictors, and modified Poisson regression applies it to a binary outcome. This article works each step by hand, then gives the Stata commands, the limits of robust standard errors and the wording for Methods.

---

🔍 Click to enlarge Click anywhere to close

Same coefficient, different standard error

A research team studies whether a treatment is associated with the number of times patients are admitted to hospital. The outcome is a count, so the team fits a Poisson regression, the standard regression for counts, which models the log of the mean count. It runs the model twice: once in the conventional way, and once with a robust standard error. After the analysis, the results are as follows.

One model, two standard errors

Variance and standard error: two levels of variability

Before a robust standard error makes sense, variability has to be split into two levels.

Level 1: the variance of the outcome. It tells how much the outcome differs between patients. For example, some patients have only 1 hospital admission, while others have 9.

Level 2: the variance of the estimator. It tells how much the estimated coefficient would change if data were collected again and again from the same population and the coefficient were estimated anew each time. The standard error is the square root of this variance:

$$\mathrm{SE}(\hat\beta) = \sqrt{\mathrm{Var}(\hat\beta)}$$

Imagine a large population of patients, and the same study repeated 1,000 times, each time on a newly sampled group of patients. If the regression coefficients from the different studies are close to one another, the estimator is precise and the standard error is low. If the coefficients swing widely between studies, the standard error is high. This is what sampling uncertainty means.

The key point: the variance of the outcome and the standard error of a coefficient are not the same thing, but how much the outcome varies does affect how uncertain the coefficient estimate is.

Two groups of patients with the same mean

Hand example. Suppose two studies each record the number of hospital admissions for 4 patients.

Same mean, different spread

Both datasets have the same mean.

$$\text{Dataset A: } (4 + 5 + 5 + 6)/4 = 5$$$$\text{Dataset B: } (1 + 1 + 9 + 9)/4 = 5$$

But look at how the data are spread. In dataset A the patients have similar numbers of admissions. In dataset B, by contrast, some patients have very few admissions and some have very many. So although the means are equal, the variance of dataset B is much higher.

Calculating the observed variance

To match the HC0 form of the sandwich estimator used later, use the empirical variance, which divides by $n$, rather than the sample variance, which divides by $n - 1$. The formula is:

$$\text{Variance} = \sum_i (Y_i - \text{Mean})^2 / n$$

For dataset A, the squared distances from the mean are 1, 0, 0 and 1:

$$\text{Variance} = (1 + 0 + 0 + 1)/4 = 0.5$$

For dataset B, every count is 4 away from the mean, so each squared distance is 16:

$$\text{Variance} = (16 + 16 + 16 + 16)/4 = 16$$

How Poisson regression sees the two datasets

Start with a Poisson regression that has only an intercept. The model is:

$$\log(\mu_i) = \beta_0$$

where $\mu_i = E(Y_i)$, the expected count for patient $i$, and log is the natural logarithm. The model has no predictors other than the intercept. Fitting it gives:

$$\hat\beta_0 = \log(\text{Mean})$$

Both datasets have a mean of 5, so:

$$\hat\beta_0 = \log(5) = 1.609$$

So the regression coefficients of dataset A and dataset B are the same.

But the conventional Poisson model also assumes something about how the data are distributed:

$$\mathrm{Var}(Y_i) = E(Y_i) = \mu_i$$

That is, when the model expects a mean of 5, it also uses a variance of 5. Compare that with what the data show.

Underdispersion and overdispersion

Compared with the Poisson variance assumption, dataset A shows underdispersion, less spread than the model assumes, and dataset B shows overdispersion, more spread than the model assumes. Even so, the conventional Poisson model computes its model-based standard error under the same assumption, variance equals mean, for both datasets.

The conventional Poisson standard error in numbers

The conventional standard error is the model-based standard error: the one implied by the variance the model assumes. For a Poisson regression with only an intercept, the model-based standard error of the estimated log mean is:

$$\mathrm{SE}_{\text{Poisson}} = \sqrt{1/(n \times \hat\mu)}$$

Here $n = 4$ patients and $\hat\mu = 5$.

1. Substitute the values

   \[ \mathrm{SE}_{\text{Poisson}} = \sqrt{1/(4 \times 5)} \]

   Four patients and a fitted mean of 5.

2. Simplify inside the root

   \[ \mathrm{SE}_{\text{Poisson}} = \sqrt{0.05} \]

   One divided by 20 is 0.05.

3. Take the square root

   \[ \mathrm{SE}_{\text{Poisson}} = 0.224 \]

   This is the standard error of the estimated log mean.

Result: The conventional standard error is 0.224 for both dataset A and dataset B, because their means are equal.

Why do both datasets get the same standard error?

This raises an important question. The counts in dataset A are close together, while those in dataset B are clearly more spread out. Why do the two datasets have the same model-based standard error?

The answer is that the conventional Poisson model does not use the observed squared residuals to estimate the standard error, as the robust sandwich estimator does. It uses the variance structure set by the Poisson model. A residual is a patient's observed count minus the fitted mean, $Y_i - \hat\mu$.

How a robust standard error fixes this

A robust standard error still uses Poisson regression to estimate the coefficient, exactly as before. What it changes is how the variance of the coefficient is calculated: it relies on the observed score contributions. A score contribution is one patient's share of the equation the software solves to find the coefficient; in a model with only an intercept, it is the patient's residual, $Y_i - \hat\mu$.

For an intercept-only Poisson model, the robust standard error in its basic form, labelled HC0, can be written as:

$$\mathrm{SE}_{\text{Robust}} = \sqrt{\sum_i (Y_i - \hat\mu)^2} \,/\, (n \times \hat\mu)$$

Now compute it for both datasets.

The robust (HC0) standard error for both datasets

Hand example. Dataset A is 4, 5, 5 and 6; dataset B is 1, 1, 9 and 9. Each has four patients and a fitted mean of 5.

1. Dataset A: sum of the squared residuals

   \[ 1 + 0 + 0 + 1 = 2 \]

   The residuals are -1, 0, 0 and +1.

2. Dataset A: robust standard error

   \[ \mathrm{SE}_{\text{Robust}} = \sqrt{2}/(4 \times 5) = 0.071 \]

   This is below the conventional 0.224.

3. Dataset B: sum of the squared residuals

   \[ 16 + 16 + 16 + 16 = 64 \]

   The residuals are -4, -4, +4 and +4.

4. Dataset B: robust standard error

   \[ \mathrm{SE}_{\text{Robust}} = \sqrt{64}/(4 \times 5) = 0.400 \]

   This is above the conventional 0.224.

Result: The robust (HC0) standard error is 0.071 for dataset A and 0.400 for dataset B.

The heart of a robust standard error

What is the sandwich variance estimator?

Now the formula of the sandwich estimator can be read with understanding [1, 2]. In a regression with several predictors, the variance of the coefficients is no longer a single number but a variance-covariance matrix. The sandwich formula is:

$$\hat V_{\text{Robust}}(\hat\beta) = A^{-1} B A^{-1}$$

It has two main parts, the bread $A$ and the meat $B$. Here the letters name the parts of the formula, not datasets A and B.

The bread: A

For Poisson regression:

$$A = \sum_i \hat\mu_i X_i X_i^\top$$

Here $X_i$ is the column of patient $i$'s predictor values, with a 1 for the intercept, and $X_i^\top$ is the same values laid out as a row. The bread describes how sensitive the estimating equations, the equations the software solves to find the coefficients, are to a change in the coefficients. Put simply, it tells how much a small change in a coefficient affects the model.

The meat: B

Start from the score contribution of observation $i$:

$$s_i = X_i (Y_i - \hat\mu_i)$$

Then compute:

$$B = \sum_i s_i s_i^\top$$

or, equivalently:

$$B = \sum_i X_i X_i^\top (Y_i - \hat\mu_i)^2$$

The meat reflects how much the score contributions vary in the observed data.

Putting the sandwich together

Combine the bread and the meat:

$$\hat V_{\text{Robust}} = A^{-1} B A^{-1}$$

Then take the square root of each value on the diagonal of the matrix to get the standard error of each coefficient:

$$\mathrm{SE}(\hat\beta_j) = \sqrt{[\hat V_{\text{Robust}}]_{jj}}$$

Back to dataset B

Hand example. With only an intercept, $X_i = 1$ for every patient, so $A$ and $B$ are single numbers that can be computed directly.

1. The bread

   \[ A = 4 \times 5 = 20 \]

   The number of patients times the fitted mean.

2. The meat

   \[ B = 16 + 16 + 16 + 16 = 64 \]

   The sum of the squared residuals.

3. The sandwich variance

   \[ \hat V_{\text{Robust}} = 64/20^2 = 0.16 \]

   The meat divided by the bread squared.

4. The robust standard error

   \[ \mathrm{SE}_{\text{Robust}} = \sqrt{0.16} = 0.40 \]

   The same value as the formula in the previous section.

Result: The sandwich gives dataset B a robust standard error of 0.40, the same HC0 value as before.

Why it is called a sandwich

$A^{-1}$ is the bread on both sides, and $B$ is the meat in the middle. In a regression with several predictors the method is the same, but the single numbers become matrix operations.

Importantly, the robust standard errors of different coefficients need not change by the same ratio, because their score contributions and the structure of their predictors differ.

How a robust standard error changes the confidence interval

Return to the regression coefficient of 0.40 from the opening table. On the log scale, using the large-sample normal approximation, the 95% confidence interval is:

$$\text{95\% CI} = \hat\beta \pm 1.96 \times \mathrm{SE}(\hat\beta)$$

1. Conventional Poisson, lower limit

   \[ 0.40 - 1.96 \times 0.10 = 0.204 \]

   Here the standard error is 0.10.

2. Conventional Poisson, upper limit

   \[ 0.40 + 1.96 \times 0.10 = 0.596 \]

   The interval runs from 0.204 to 0.596.

3. Robust Poisson, lower limit

   \[ 0.40 - 1.96 \times 0.19 = 0.028 \]

   Here the standard error is 0.19.

4. Robust Poisson, upper limit

   \[ 0.40 + 1.96 \times 0.19 = 0.772 \]

   The interval runs from 0.028 to 0.772.

5. Exponentiate the coefficient

   \[ \exp(0.40) \approx 1.49 \]

   The same for both methods.

Result: Both methods give the same point estimate, 0.40, or 1.49 after exponentiating; the larger robust standard error makes the confidence interval wider.

Linking to modified Poisson regression for a binary outcome

In clinical epidemiology the outcome is often binary, such as neonatal morbidity within 28 days: $Y = 1$ when the outcome occurs and $Y = 0$ when it does not. For a binary outcome with probability $p$ of the event, the expected value is:

$$E(Y) = p$$

But the variance under the Bernoulli distribution, the Bernoulli variance, is:

$$\mathrm{Var}(Y) = p(1 - p)$$

while the Poisson variance assumption is:

$$\mathrm{Var}(Y) = p$$

Suppose the risk is 0.40. The Poisson model uses a variance of 0.40, but the Bernoulli variance is:

$$0.40 \times (1 - 0.40) = 0.24$$

So the Poisson variance assumption does not match the Bernoulli variance. Importantly, this mismatch is not overdispersion, as in the count example earlier: the Bernoulli variance is below the Poisson variance whenever $0 < p < 1$.

The model and what it estimates

The model is:

$$\log[P(Y = 1 \mid X)] = \beta_0 + \beta_1 X$$

When $X$ is a binary exposure and the model is appropriate:

$$\mathrm{RR} = \exp(\beta_1)$$

So modified Poisson regression does not turn the outcome into real count data. It applies the Poisson estimating framework to a binary outcome.

It is one alternative to log-binomial regression, which models the log of the risk directly, especially when the log-binomial model has convergence problems. However, a robust standard error does not fix predicted risks that may exceed 1, nor a mean model that is specified wrongly. The article on risk ratios for a common outcome compares the two models in more depth.

How Stata computes a robust standard error

The commands below are written for the admissions example from the opening: a count of hospital admissions, a treatment indicator and age. Stata sets how the variance is estimated through its vce() option [4].

Conventional Poisson, then Poisson with a robust standard error

A. Conventional Poisson:

poisson admission_count i.treatment c.age

B. Poisson with a robust standard error:

poisson admission_count i.treatment c.age, vce(robust)

With the same data and the same model specification, the two commands give the same regression coefficients. Their standard errors, confidence intervals and p-values may differ. The option vce(robust) asks Stata to use the robust variance-covariance estimator.

Modified Poisson for a binary outcome

C. To study the association between inadequate antenatal care (ANC) and neonatal morbidity:

poisson morbidity28 i.inadequate_anc i.age_group i.low_ses i.coverage_gap, irr vce(robust)

If the outcome is binary, everyone is followed for the same period, and the log-risk model is appropriate, the exponentiated coefficient can be read as an adjusted risk ratio. Stata heads the column IRR because of the irr option, short for incidence rate ratio, but in this setting the effect measure is the risk ratio, not an incidence rate ratio.

By contrast, if the outcome is a number of events per person-time and the exposure time enters the model, the exponentiated coefficient means a rate ratio.

So the name of the effect measure must come from the outcome, the time denominator and the estimand, not just from the name of a Stata option. The estimand is the precise quantity the study aims to estimate.

Cluster-robust standard errors

D. When patients have repeated measurements:

poisson admission_count i.treatment c.age, vce(cluster patient_id)

This command computes a cluster-robust standard error: it adds up the score contributions within the same patient before building the meat of the sandwich estimator. Large-sample (asymptotic) inference needs enough independent clusters, and when the clusters are few, a method that corrects the small-sample bias should be considered. The same idea extends modified Poisson regression to correlated binary data [5]. Part 2 of this series works through when to use which vce() option.

What a robust standard error can and cannot do

Overdispersed counts: what is the analysis for?

When a count outcome is overdispersed in particular, the researcher has to consider what the analysis is for. If the interest is the conditional mean association and the mean model is correct, Poisson regression with a robust standard error may be enough for inference. If the interest is the count distribution, a predictive distribution or a description of the variance structure, consider negative binomial regression, a count model that lets the variance exceed the mean, or another suitable model.

A robust standard error does not make the variance of the data smaller, and it does not automatically turn the regression model into a correct one [6].

How to write it in the Methods of a manuscript

For a binary outcome in a cohort study, it can be written as:

"We estimated adjusted risk ratios using modified Poisson regression with a log link and robust sandwich variance estimation to obtain 95% confidence intervals."

State the adjustment variables, the handling of missing data, the follow-up and any clustering, as the study design requires.

For a count outcome, it can be written as:

"We fitted Poisson regression models with robust sandwich standard errors to estimate associations with the expected event count."

If a person-time offset is included and the interest is the incidence rate ratio, also state that the log of person-time was used as the offset. An offset is a term entered with its coefficient fixed at 1, so that the model describes events per unit of time.

Conclusion: point estimation and variance estimation

The most important thing about robust standard errors is not memorising the sandwich formula. It is understanding that a statistical model has two parts that must be kept apart.

Part one: point estimation. The model uses the data to estimate the regression coefficient, which gives the size and direction of the association.

Part two: variance estimation. A statistical method assesses how uncertain that estimated coefficient is.

Conventional Poisson regression uses the variance assumption of the Poisson model. Robust Poisson regression uses the sandwich estimator to estimate the uncertainty from the observed score contributions, while keeping the same mean model.

So a robust standard error is not a technique that makes the data better. It is a tool that makes inference about a coefficient more appropriate when the variance structure of the model does not match the data, under the assumptions and conditions it needs.

Common misreadings and their fixes

• "A robust standard error means a new model."

  The coefficient in the opening table stays at 0.40 with either standard error, and both hand datasets keep the coefficient 1.609.

  Fix: A robust standard error keeps the same model and the same coefficient; it changes only how the variance-covariance matrix is estimated.

• "A robust standard error is always larger."

  In dataset A the robust standard error is 0.071, below the conventional 0.224.

  Fix: The robust standard error follows the observed spread: smaller when the data vary less than the model assumes, larger when they vary more.

• "Poisson regression of a binary outcome suffers from overdispersion."

  At a risk of 0.40 the Bernoulli variance is 0.24, below the Poisson variance of 0.40.

  Fix: The mismatch runs the other way: the Bernoulli variance is below the Poisson variance whenever the risk lies between 0 and 1, and the robust standard error takes its place in the inference.

• "Stata prints IRR, so the result is an incidence rate ratio."

  The irr option only sets the column label. With a binary outcome and the same follow-up for everyone, the exponentiated coefficient is a risk ratio.

  Fix: Name the effect measure from the outcome, the time denominator and the estimand, not from the Stata option.

• "A robust standard error makes the model correct."

  It does not fix a wrongly specified mean model, confounding, selection bias, an unspecified cluster structure or a very small sample, and it does not turn a Poisson model into a negative binomial one.

  Fix: Use a robust standard error for the uncertainty of the coefficient; deal with the mean model, confounding and the design separately.

Glossary

Poisson regression

A regression for counts that models the log of the mean count and assumes that the variance equals the mean.

variance-covariance matrix (เมทริกซ์ความแปรปรวนและความแปรปรวนร่วม)

The matrix holding each coefficient's estimated variance on its diagonal and the covariances between coefficients off it; standard errors are the square roots of the diagonal.

sampling uncertainty

How much an estimate would vary between repeated studies on new samples from the same population; the standard error measures it.

empirical variance

The average squared distance of the observations from their mean, dividing by n rather than n - 1.

underdispersion (การกระจายน้อยกว่าที่แบบจำลองกำหนด)

Less spread in a count outcome than the Poisson model assumes: the variance is below the mean.

overdispersion (การกระจายเกินกว่าที่แบบจำลองกำหนด)

More spread in a count outcome than the Poisson model assumes: the variance exceeds the mean.

model-based standard error (standard error แบบอิงแบบจำลอง)

The standard error implied by the variance the model assumes, such as variance equal to the mean in Poisson regression.

score contribution (score contribution ของแต่ละหน่วยสังเกต)

One observation's share of the estimating equations: its predictor values times its residual; with only an intercept, the residual itself.

HC0

The basic robust (sandwich) standard error, without the finite-sample correction that some software applies.

sandwich estimator

A variance estimate built as inverse bread, meat, inverse bread; it keeps the coefficient and reads the spread from the observed score contributions.

bread (ส่วน bread ของ sandwich estimator)

The outer factor of the sandwich, A: how sensitive the estimating equations are to the coefficients; it appears inverted on both sides of the meat.

meat (ส่วน meat ของ sandwich estimator)

The middle factor of the sandwich, B: the sum of the squared score contributions, read from the data rather than from the variance the model assumes.

Bernoulli variance (ความแปรปรวนแบบ Bernoulli)

The variance p(1 - p) of a yes-or-no outcome with probability p; it is below the Poisson value p whenever p lies between 0 and 1.

modified Poisson regression (Poisson ดัดแปลง)

A Poisson log-link mean model for a binary outcome, with a robust sandwich standard error, used to estimate a risk ratio.

log-binomial regression (แบบจำลอง log-binomial)

A regression for a yes-or-no outcome that models the log of the risk directly; it can fail to converge.

rate ratio (อัตราส่วนอัตราอุบัติการณ์)

The ratio of events per person-time between two groups; Stata labels it IRR, short for incidence rate ratio.

estimand (ปริมาณเป้าหมายของการประมาณ)

The precise quantity a study aims to estimate, including the population, contrast and summary measure.

cluster-robust standard error (cluster-robust standard error)

A sandwich standard error whose meat adds up the score contributions within each cluster, such as one patient, before squaring; it needs enough independent clusters.

negative binomial regression (negative binomial regression)

A count model with the Poisson mean structure plus a dispersion parameter, so the variance can exceed the mean.

offset

A term entered into the model with its coefficient fixed at 1, such as the log of person-time, so the model describes events per unit of time.

References

1. Huber PJ. The behavior of maximum likelihood estimates under nonstandard conditions. In: Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1. University of California Press; 1967. p. 221-233. https://projecteuclid.org/euclid.bsmsp/1200512988
2. White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838. doi:10.2307/1912934 https://doi.org/10.2307/1912934
3. Zou G. A modified Poisson regression approach to prospective studies with binary data. Am J Epidemiol. 2004;159(7):702-706. doi:10.1093/aje/kwh090 https://doi.org/10.1093/aje/kwh090
4. StataCorp. Stata Base Reference Manual: [R] vce_option. Stata documentation at stata.com (PDF). Accessed 2026-10-11. https://www.stata.com/manuals/rvce_option.pdf
5. Zou GY, Donner A. Extension of the modified Poisson regression model to prospective studies with correlated binary data. Stat Methods Med Res. 2013;22(6):661-670. doi:10.1177/0962280211427759 https://doi.org/10.1177/0962280211427759
6. Freedman DA. On the so-called 'Huber sandwich estimator' and 'robust standard errors'. Am Stat. 2006;60(4):299-302. doi:10.1198/000313006X152207 https://doi.org/10.1198/000313006X152207

Key takeaways

• Estimating a coefficient and estimating its uncertainty are separate processes: a robust standard error keeps the model and the coefficient and changes only the variance-covariance estimate.
• With the same mean, datasets A and B share the coefficient 1.609 and the conventional standard error 0.224, but their robust (HC0) standard errors are 0.071 and 0.400, because their observed spread differs.
• The sandwich is A⁻¹BA⁻¹: the bread A comes from the model, the meat B from the observed score contributions, and each standard error is the square root of a diagonal value.
• Modified Poisson regression uses a Poisson log-link model for the risk ratio and a robust sandwich variance for inference; name the effect measure from the outcome, the time denominator and the estimand, not from the Stata option.
• Robust standard errors change how uncertainty is estimated, not what coefficient is estimated; they do not fix a wrong mean model, confounding or selection bias.

Related in the wiki: [[heteroskedasticity-robust-standard-errors]] [[risk-regression-models-epidemiology]] [[choosing-vce-robust-cluster-standard-errors]]