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

Published by UniqcretEN12 October 2026· 25 min readอ่านภาษาไทย
Clinical Epidemiology ResearchMethodology and Research DesignUniqcret doctor knowledges
Same Coefficient, Different Standard Error: The Sandwich Estimator Worked by Hand
On this page

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

Abstract

A Poisson regression fitted twice, plain and with a robust standard error, prints one coefficient beside two standard errors. The coefficient comes from the mean model; the standard error comes from a separate variance calculation. Two hand datasets of four patients, too few for inference, share a mean of 5 admissions. An intercept-only Poisson model gives both the coefficient ln 5 = 1.609 and the model-based standard error 0.224. The sandwich estimator reads the spread from the data instead. Its bread, J = 20, is the same in both datasets. Its meat, K, is 2 in one and 64 in the other. The HC0 sandwich standard errors are 0.071 and 0.400. The sandwich shrinks the standard error where counts vary less than assumed, and widens it where they vary more. Stata's vce(robust) also multiplies the variance by n/(n - 1), so it prints larger values. This article works each step by hand, checks it in Stata and R, and extends it to a binary outcome.

Visual summary. Hand example.

One model, two standard errors

A hospital pharmacy team audits admissions across a regional network to plan its stock. Its file holds one day's admission count for each of 100 hospitals. The counts swing from 0 to 10 around a mean of 5 admissions per day. The hospitals are fictional and the data are simulated.

The analyst fits a Poisson regression, the standard model for counts, which models the log of the mean count. With no covariates, its one coefficient is the log of the mean. The analyst runs it twice in Stata: once plain, and once with the option vce(robust). That option asks Stata for a robust (sandwich) estimate of the variance.

The coefficient is identical in both runs, while the standard error almost doubles, as the table below shows. The plain run reports the model-based standard error, which trusts the variance the Poisson model assumes. The team wants to know how one model can be sure and unsure of the same number at once. This article works it out by hand on four patients, then returns to the 100 hospitals.

The same coefficient with two standard errors

Simulated data, 100 hospitals, one day each. The coefficient and the mean do not move between the runs; the standard error and both intervals do.
QuantityModel-based runRun with vce(robust)
Coefficient, ln(mean admissions per day)1.60941.6094
Standard error of the coefficient0.04470.0853
95% confidence interval (CI) of the coefficient1.5218 to 1.69711.4423 to 1.7766
Mean admissions per day, exp(coefficient)55
95% CI of the mean (admissions per day)4.58 to 5.464.23 to 5.91

Two kinds of spread

Two different variances meet in this problem. The variance of the outcome is how much admission counts differ between patients. The variance of the estimate is how much a coefficient would differ between repeated studies.

Picture the same study run 1,000 times, each time on a fresh sample from the same population. The standard error is, in effect, the standard deviation of those 1,000 coefficients: the square root of the second variance. A single study never sees the other runs, so it must estimate that spread from its own data.

The two variances are linked. The more the counts vary between patients, the more an estimate built from them varies between studies. The two standard errors in this article differ only in where that outcome variance comes from. One takes it from the model's assumption, the other from the data.

Hand example: two datasets with the same mean

Hand example. Two datasets each record the hospital admissions of four patients over the same period. Dataset A has counts of 4, 5, 5 and 6; dataset B has 1, 1, 9 and 9. Each dataset adds up to 20 admissions, so both means are 5.

Dataset A is tightly packed, while dataset B splits into two very different pairs. The empirical variance measures that spread: the average squared distance from the mean, dividing by $n$, the number of patients. Dividing by $n$ rather than $n - 1$ matches the sandwich formula used below. In dataset A the squared distances from the mean are 1, 0, 0 and 1.

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

So dataset A's empirical variance is 0.5. In dataset B every count sits 4 admissions from the mean, so each squared distance is 16.

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

So dataset B's empirical variance is 16, far above dataset A's 0.5.

What the Poisson model assumes

An intercept-only Poisson regression, one with no covariates, fits the log of the mean count:

$$\hat\beta_0 = \ln \bar y = \ln 5 = 1.609$$

Here $\bar y$ is the mean count and $\hat\beta_0$ is the estimated intercept; a hat marks an estimate. Both datasets have a mean of 5, so both get the same coefficient. The fitted mean, $\hat\mu = 5$, is the same for every patient in a dataset.

The Poisson model also assumes that the variance of a count equals its mean, here 5. Dataset A, with an empirical variance of 0.5, shows underdispersion: less spread than the model assumes. Dataset B, at 16, shows overdispersion: more spread than the model assumes.

The model-based standard error

The model-based standard error is the one implied by the variance the model assumes. For an intercept-only Poisson model with $n$ patients and fitted mean $\hat\mu$, it is:

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

Here $n\hat\mu = 4 \times 5 = 20$ is the total expected count. The standard error is $\sqrt{1/20} = \sqrt{0.05} = 0.224$ on the log scale, for both datasets. The formula sees only $n$ and $\hat\mu$, so it cannot tell dataset A from dataset B.

The sandwich: bread and meat

The sandwich estimator keeps the coefficient but recomputes its variance from the data [1, 2]. Its standard error is often called the robust standard error. Its basic form, with no small-sample factor, is labelled HC0. For an intercept-only Poisson model, each piece of the sandwich is a single number.

The score contribution of patient $i$ is $s_i = y_i - \hat\mu$, the patient's count $y_i$ minus the fitted mean. With covariates, it is that residual times the patient's covariate values. The meat, $K$, adds these up, squared:

$$K = \sum_i s_i^2 = \sum_i (y_i - \hat\mu)^2$$

So the meat is read from the data: each patient adds their squared distance from the fitted mean. The bread, $J$, measures how sharply the estimating equation changes as the coefficient moves. That equation is the one the software solves to find the coefficient. For this model the bread is:

$$J = n\hat\mu$$

That is the total expected count, the product inside the model-based formula. The sandwich variance puts the inverse of the bread on both sides of the meat:

$$\widehat{\mathrm{Var}}(\hat\beta_0) = J^{-1} K J^{-1} = K/J^2$$

Its square root is the HC0 standard error. The model-based variance is the inverse bread alone, $1/(n\hat\mu)$. If the counts vary as the Poisson model assumes, $K$ lands close to $J$ and the two standard errors nearly agree. With covariates, $J$ and $K$ become matrices, as the article on robust standard errors shows.

Hand example: the bread and the meat as numbers

Hand example. Dataset A has counts of 4, 5, 5 and 6, and dataset B has 1, 1, 9 and 9. Each has four patients and a fitted mean of $\hat\mu = 5$.

  1. The bread, the same for both datasets

    \[ J = n\hat\mu = 4 \times 5 = 20 \]

    It depends only on the number of patients and the fitted mean.

  2. The meat for dataset A

    \[ K = (-1)^2 + 0^2 + 0^2 + 1^2 = 2 \]

    The score contributions are -1, 0, 0 and +1.

  3. The meat for dataset B

    \[ K = 4 \times 16 = 64 \]

    Each score contribution is -4 or +4, so each squares to 16.

  4. The bread squared

    \[ J^2 = 20^2 = 400 \]

    This is the same for both datasets.

  5. Sandwich variance for dataset A

    \[ K/J^2 = 2/400 = 0.005 \]

    A small meat gives a small variance.

  6. Sandwich variance for dataset B

    \[ K/J^2 = 64/400 = 0.16 \]

    A large meat gives a large variance.

  7. HC0 standard error for dataset A

    \[ \sqrt{0.005} = 0.071 \]

    This is below the model-based standard error of 0.224.

  8. HC0 standard error for dataset B

    \[ \sqrt{0.16} = 0.400 \]

    This is above the model-based standard error of 0.224.

Result: Both datasets keep the coefficient 1.609 and the model-based standard error 0.224. Their HC0 standard errors, 0.071 and 0.400, differ because only the meat sees how the counts spread.

A shortcut that holds only for this model

For the intercept-only Poisson model, the steps collapse into one line:

$$\mathrm{SE}_{\text{HC0}} = \sqrt{K}/J = \sqrt{K}/(n\hat\mu)$$

For dataset B this gives $\sqrt{64}/20 = 8/20 = 0.400$. The shortcut holds for the intercept-only model only, because only there is the bread the single number $n\hat\mu$. With covariates, each standard error comes from the diagonal of the matrix $J^{-1} K J^{-1}$.

Hand example. Edit the four admission counts, or load dataset A or B with the buttons. The panel recomputes the mean, ln(mean), the bread $J$, the meat $K$, both standard errors and their ratio. It starts at dataset B. Counts that vary about as much as their mean bring the two standard errors together. Four counts show the arithmetic only. Four equal counts give an HC0 standard error of 0, a sign of too few patients, not of perfect precision.

Which way the sandwich moves

In this intercept-only model, divided by $n$, the meat is the empirical variance. The bread, divided by $n$, is the variance the model assumes. So here the HC0 standard error exceeds the model-based one exactly when $K$ exceeds $J$.

Dataset A has $K = 2$ against $J = 20$, so the sandwich shrinks its standard error. Dataset B has $K = 64$ against $J = 20$, so the sandwich widens it. With covariates, $J$ and $K$ are matrices, and no single comparison of the two decides the direction. It depends on where the large residuals sit, as the article on robust standard errors sets out.

Four patients are enough to show the arithmetic but far too few to trust either standard error. The sandwich is a large-sample method, and with few observations its HC0 value tends to be too small. That article also compares versions that correct for this.

The same arithmetic in Stata and R

The series' simulation script types both datasets in as fixed counts and repeats each step with matrix commands. In Stata, matrix accum builds $J$ with each patient weighted by $\hat\mu$, and $K$ with each weighted by the squared residual. The panes show dataset B. For dataset A, the same lines print $J = 20$, $K = 2$ and a sandwich variance of 0.005.

Stata: the model-based fit, then the bread and the meat

Stata code g1_sim.do (lines 57-70 of 257)
    * model-based standard error: the Poisson model assumes variance = mean
    poisson admissions if dataset == `d', nolog
    canon `k'.b0 4 _b[_cons]
    canon `k'.se_model 4 _se[_cons]
    * bread and meat by hand (one column, the constant)
    predict double mu if dataset == `d', n
    generate double r2 = (admissions - mu)^2 if dataset == `d'
    matrix accum J = one if dataset == `d' [iweight = mu], noconstant
    matrix accum K = one if dataset == `d' [iweight = r2], noconstant
    matrix Ji = invsym(J)
    matrix V = Ji*K*Ji
    matrix list J
    matrix list K
    matrix list V
Output of the run g1_sim.log
Poisson regression                                      Number of obs =      4
                                                        LR chi2(0)    =   0.00
                                                        Prob > chi2   =      .
Log likelihood = -13.414897                             Pseudo R2     = 0.0000

------------------------------------------------------------------------------
  admissions | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .2236068     7.20   0.000     1.171177    2.047699
------------------------------------------------------------------------------
CANON g1.b.b0 1.6094
CANON g1.b.se_model 0.2236
(4 missing values generated)
(4 missing values generated)
(obs=20)
(obs=64)

symmetric J[1,1]
     one
one   20

symmetric K[1,1]
     one
one   64

symmetric V[1,1]
     one
one  .16
Hand example data, an excerpt from the series' simulation script and its output. The file name and line numbers only record where the excerpt sits. The lines run once per dataset, and the output shown is dataset B's. It gives the model-based standard error, 0.2236, then $J = 20$, $K = 64$ and the matrix V, the sandwich variance, 0.16. In the code, Ji is the inverse of $J$. (obs=20) and (obs=64) are the weighted totals $J$ and $K$. The two missing-value lines appear because predict and generate fill only the current dataset's rows. Lines beginning CANON come from a small helper in the script that echoes each result; read the number at the end.

R: the same fit, the bread and the meat, and vcovHC()

R code g1_sim_r.R (lines 55-74 of 205)
  # model-based standard error: the Poisson model assumes variance = mean
  fit <- glm(admissions ~ 1, family = poisson, data = s)
  print(summary(fit)$coefficients)
  canon(paste0(k, ".b0"), coef(fit)[[1]])
  canon(paste0(k, ".se_model"), sqrt(vcov(fit)[1, 1]))
  # bread and meat by hand (one column, the constant)
  bm <- bread_meat(matrix(1, n, 1), y, fitted(fit))
  print(bm)
  canon(paste0(k, ".bread"), bm$J[1, 1])
  canon(paste0(k, ".meat"), bm$K[1, 1])
  canon(paste0(k, ".v_model"), bm$Ji[1, 1], 6)
  canon(paste0(k, ".v_hc0"), bm$V[1, 1], 6)
  canon(paste0(k, ".se_hc0"), sqrt(bm$V[1, 1]))
  # Stata's robust standard error: the sandwich package's HC0 times n/(n - 1)
  canon(paste0(k, ".se_robust"), sqrt(sandwich(fit)[1, 1] * n / (n - 1)))
  canon(paste0(k, ".robust_factor"), n / (n - 1))
  # the same HC0 variance from sandwich::vcovHC() (for a GLM it equals sandwich(fit) and the hand value above)
  V_hc0 <- sandwich::vcovHC(fit, type = "HC0")
  print(V_hc0)
  stopifnot(isTRUE(all.equal(V_hc0[1, 1], bm$V[1, 1])))
Output of the run g1_sim_r.log
            Estimate Std. Error  z value     Pr(>|z|)
(Intercept) 1.609438  0.2236067 7.197629 6.126851e-13
CANON g1.b.b0 1.6094
CANON g1.b.se_model 0.2236
$J
     [,1]
[1,]   20

$K
     [,1]
[1,]   64

$Ji
     [,1]
[1,] 0.05

$V
     [,1]
[1,] 0.16

CANON g1.b.bread 20.0000
CANON g1.b.meat 64.0000
CANON g1.b.v_model 0.050000
CANON g1.b.v_hc0 0.160000
CANON g1.b.se_hc0 0.4000
CANON g1.b.se_robust 0.4619
CANON g1.b.robust_factor 1.3333
            (Intercept)
(Intercept)        0.16
Hand example data, an excerpt from the R version of the series' simulation script and its output for dataset B. The lines run once per dataset inside a loop; s holds one dataset's four rows. The function bread_meat(), defined earlier in the script, returns $J$, $K$, the inverse bread Ji and the sandwich variance V. Then sandwich::vcovHC(fit, type = "HC0") prints the same HC0 variance, 0.16. Leave out type = "HC0" and vcovHC() returns a different small-sample version. Lines beginning CANON come from a small helper that echoes each result. For the robust line, the script multiplies HC0 by $n/(n - 1)$ to match Stata, as the next section explains.

Why Stata prints a larger robust standard error

Fit dataset B with vce(robust) and Stata prints a standard error of 0.4619, not 0.400. The option asks Stata for a robust (sandwich) variance [3]. For commands that fit by maximum likelihood, such as poisson or logistic, Stata multiplies that variance by $n/(n - 1)$. The values worked by hand above are therefore labelled HC0, the sandwich without that factor.

With four patients the factor is large:

$$n/(n - 1) = 4/3 = 1.333$$

The variance grows by a factor of 1.333. The standard error grows by the square root of that factor:

$$\sqrt{4/3} = 1.155$$

So the standard error grows by a factor of 1.155. Dataset A becomes $0.071 \times 1.155 = 0.082$, and dataset B becomes $0.400 \times 1.155 = 0.462$. These are Stata's printed values, rounded to three decimals, so the run confirms the factor. With 100 hospitals the factor is barely above 1.

Stata: the same fit with vce(robust)

Stata code g1_sim.do (lines 76-77 of 257)
    * Stata's robust standard error: the same sandwich times n/(n - 1)
    poisson admissions if dataset == `d', vce(robust) nolog
Output of the run g1_sim.log
Poisson regression                                      Number of obs =      4
                                                        Wald chi2(0)  =      .
                                                        Prob > chi2   =      .
Log pseudolikelihood = -13.414897                       Pseudo R2     = 0.0000

------------------------------------------------------------------------------
             |               Robust
  admissions | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .4618802     3.48   0.000     .7041693    2.514706
------------------------------------------------------------------------------
Hand example data, an excerpt from the series' simulation script and its output for dataset B. The coefficient is unchanged. The robust standard error, 0.4619, is the HC0 value 0.400 times the square root of 4/3.

The interval: same centre, different width

Return to the 100 hospitals of the opening table. Their counts vary from 0 to 10, with an empirical variance of 18 against a mean of 5. The meat, $K$, is 1,800 against a bread, $J$, of 500, so the sandwich widens the standard error. The HC0 standard error is $\sqrt{1800}/500 = 0.0849$.

Stata's vce(robust) adds its factor $n/(n - 1)$ and prints 0.0853, 1.91 times the model-based 0.0447. Both 95% intervals are centred on the same coefficient, 1.6094, and the robust one is 1.91 times as wide. On the log scale it runs from 1.4423 to 1.7766, against 1.5218 to 1.6971. Exponentiating the limits gives limits for the mean count, which stays at 5 admissions per day.

A steady group of 100 other hospitals has the same mean, but counts of only 4, 5 or 6. Its empirical variance is 0.5, so its robust interval, 1.5816 to 1.6373, is narrower than its model-based one.

One model can hold both groups, with an indicator for the volatile group. The indicator's coefficient is the log of a rate ratio, the ratio of the two mean counts per day. With equal means the rate ratio is 1.00, and only the width of its interval depends on the standard error. Its robust interval is the wider one, because the volatile group's extra spread outweighs the steady group's shortfall.

Model-based and robust intervals in the two hospital groups

Simulated data, two groups of 100 hospitals, one day each, mean 5 admissions per day in both. The estimates do not move; only the 95% intervals change.
QuantityEstimateModel-based 95% CIRobust 95% CI
Volatile group, coefficient ln(mean)1.60941.5218 to 1.69711.4423 to 1.7766
Steady group, coefficient ln(mean)1.60941.5218 to 1.69711.5816 to 1.6373
Steady group, mean admissions per day54.58 to 5.464.86 to 5.14
Volatile versus steady, log rate ratio0-0.1240 to 0.1240-0.1690 to 0.1690
Volatile versus steady, rate ratio1.000.88 to 1.130.84 to 1.18

Stata: the volatile group, model-based and robust

Stata code g1_sim.do (lines 112-129 of 257)
    poisson admissions if volatile == `g', nolog
    scalar sc_b = _b[_cons]
    scalar sc_sm = _se[_cons]
    canon `k'.b0 4 _b[_cons]
    canon `k'.se_model 4 _se[_cons]
    predict double mu if volatile == `g', n
    generate double r2 = (admissions - mu)^2 if volatile == `g'
    matrix accum J = one if volatile == `g' [iweight = mu], noconstant
    matrix accum K = one if volatile == `g' [iweight = r2], noconstant
    matrix Ji = invsym(J)
    matrix V = Ji*K*Ji
    canon `k'.bread 4 el(J, 1, 1)
    canon `k'.meat 4 el(K, 1, 1)
    canon `k'.v_model 6 el(Ji, 1, 1)
    canon `k'.v_hc0 6 el(V, 1, 1)
    canon `k'.se_hc0 4 sqrt(el(V, 1, 1))
    canon `k'.ratio_hc0_model 4 sqrt(el(V, 1, 1))/scalar(sc_sm)
    poisson admissions if volatile == `g', vce(robust) nolog
Output of the run g1_sim.log
Poisson regression                                      Number of obs =    100
                                                        LR chi2(0)    =   0.00
                                                        Prob > chi2   =      .
Log likelihood = -365.2899                              Pseudo R2     = 0.0000

------------------------------------------------------------------------------
  admissions | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .0447214    35.99   0.000     1.521786     1.69709
------------------------------------------------------------------------------
CANON g1.volatile.b0 1.6094
CANON g1.volatile.se_model 0.0447
(100 missing values generated)
(100 missing values generated)
(obs=500)
(obs=1800)
CANON g1.volatile.bread 500.0000
CANON g1.volatile.meat 1800.0000
CANON g1.volatile.v_model 0.002000
CANON g1.volatile.v_hc0 0.007200
CANON g1.volatile.se_hc0 0.0849
CANON g1.volatile.ratio_hc0_model 1.8974

Poisson regression                                      Number of obs =    100
                                                        Wald chi2(0)  =      .
                                                        Prob > chi2   =      .
Log pseudolikelihood = -365.2899                        Pseudo R2     = 0.0000

------------------------------------------------------------------------------
             |               Robust
  admissions | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .0852803    18.87   0.000     1.442292    1.776584
------------------------------------------------------------------------------
Simulated data, an excerpt from the series' simulation script and its output for the volatile group. The lines run once per group of 100 hospitals. The first table prints the model-based interval, 1.5218 to 1.6971 at four decimals. The second prints the robust interval, 1.4423 to 1.7766. Between them, (obs=500) and (obs=1800) are the bread and the meat, and lines beginning CANON come from a small helper that echoes each result.

A yes-or-no outcome: modified Poisson regression

The same arithmetic explains a common choice for binary outcomes. A yes-or-no outcome with risk $p$ has the Bernoulli variance $p(1 - p)$, while a Poisson model assumes $p$, the mean. At $p = 0.40$ the Bernoulli variance is $0.40 \times 0.60 = 0.24$, against the Poisson 0.40. For any risk between 0 and 1, $p(1 - p)$ is below $p$.

A Poisson model therefore overstates the variance of a binary outcome. Its model-based standard error tends to be too large, and the sandwich reads the smaller spread from the data. Modified Poisson regression uses exactly this pair [4]. It fits a Poisson model with a log link for the risk ratio and takes the standard error from the sandwich.

The simulated dataset holds 200 newborns, half with inadequate antenatal care (ANC). Neonatal morbidity within four weeks of birth has a risk of 0.40 after inadequate ANC and 0.20 after adequate ANC. Both Poisson fits give a risk ratio of 2.00, a coefficient of 0.6931 on the log scale. The model-based standard error is 0.2739; the sandwich gives 0.2345 as HC0 and 0.2351 from vce(robust).

Here the robust standard error is the smaller one, as the Bernoulli arithmetic predicts. A log-binomial model models the log of the risk and builds in the Bernoulli variance. It gives the same risk ratio, 2.00, with a standard error of 0.2345, the HC0 value. On the log scale, the Poisson coefficient's 95% interval narrows from 0.1564 to 1.2299 (model-based) to 0.2323 to 1.1540 (robust).

Logistic regression gives an odds ratio of 2.67 instead, because it models the log odds. Both ratios are unadjusted comparisons of the same 200 newborns. The odds ratio sits further from 1 because the outcome is common, not because of confounding. The article on risk ratios for a common outcome sets out when to choose modified Poisson and what to try first.

One binary outcome, three models

Simulated data, 200 newborns. The Poisson fits share one coefficient, and the sandwich narrows its interval. Logistic regression estimates an odds ratio, a different measure.
ModelMeasureEstimateStandard error of the coefficient95% CI
Poisson, model-basedRisk ratio2.000.27391.17 to 3.42
Poisson, HC0 sandwich by matrix commandsRisk ratio2.000.23451.26 to 3.17
Modified Poisson, vce(robust)Risk ratio2.000.23511.26 to 3.17
Logistic regressionOdds ratio2.670.32271.42 to 5.02

Stata: modified Poisson with vce(robust)

Stata code g1_sim.do (lines 240-241 of 257)
* modified Poisson: the same coefficient with Stata's robust (sandwich) standard error
poisson morbidity inadequate_anc, irr vce(robust) nolog
Output of the run g1_sim.log
. * modified Poisson: the same coefficient with Stata's robust (sandwich) standard error
. poisson morbidity inadequate_anc, irr vce(robust) nolog

Poisson regression                                      Number of obs =    200
                                                        Wald chi2(1)  =   8.69
                                                        Prob > chi2   = 0.0032
Log pseudolikelihood = -128.84039                       Pseudo R2     = 0.0257

--------------------------------------------------------------------------------
               |               Robust
     morbidity |        IRR   std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |          2   .4702186     2.95   0.003      1.26155    3.170702
         _cons |         .2   .0401004    -8.03   0.000     .1350089    .2962767
--------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.
Simulated data, an excerpt from the series' simulation script and its output. With irr, Stata prints the ratio, 2. Its standard error column refers to the ratio itself: the coefficient's robust standard error, 0.2351, times 2. The header says IRR, Stata's label for an incidence rate ratio; with a binary outcome and equal follow-up, the ratio is a risk ratio. The _cons row, 0.2, is the risk after adequate ANC, whatever Stata's note below the table calls it.

Stata: logistic regression for contrast

Stata code g1_sim.do (lines 248-249 of 257)
* logistic regression for contrast: an odds ratio, not a risk ratio
logistic morbidity inadequate_anc, nolog
Output of the run g1_sim.log
. * logistic regression for contrast: an odds ratio, not a risk ratio
. logistic morbidity inadequate_anc, nolog

Logistic regression                                     Number of obs =    200
                                                        LR chi2(1)    =   9.66
                                                        Prob > chi2   = 0.0019
Log likelihood = -117.34141                             Pseudo R2     = 0.0395

--------------------------------------------------------------------------------
     morbidity | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |   2.666666   .8606626     3.04   0.002     1.416591    5.019872
         _cons |   .2500001      .0625    -5.55   0.000      .153158    .4080755
--------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
Simulated data, an excerpt from the series' simulation script and its output. The same outcome and exposure give an odds ratio of 2.67, a different measure from the risk ratio of 2.00.

What the sandwich repairs, and what comes next

Every example in this article separates two jobs. The mean model sets the coefficient; the variance estimate sets its standard error and interval. The sandwich changes only the second job.

A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias [5].

The patients and hospitals here were independent of one another. When observations are grouped, for example patients within hospitals, the choice of vce() changes. Part 2 of this series works through that choice and how to report it.

Common misreadings and their fixes

  • "A robust standard error is always larger."

    For the binary outcome in 200 newborns, the robust standard error, 0.2351, is below the model-based 0.2739. The hand example points the same way: dataset A's HC0 standard error is 0.071 against 0.224. Four patients, though, only show the arithmetic.

    Fix: The sandwich can move the standard error either way: here it is smaller when the counts vary less than the Poisson model assumes and larger when they vary more.

  • "Robust standard errors fix a misspecified model."

    The sandwich leaves the coefficient where the mean model put it. A wrong mean model, such as one that leaves out a confounder, gives a wrong coefficient under every variance estimator.

    Fix: A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias.

  • "Stata's vce(robust) prints the HC0 sandwich."

    For a maximum likelihood command such as poisson, Stata multiplies the HC0 variance by n/(n - 1). With four patients, that raises the standard error by a factor of 1.155.

    Fix: Label hand-computed values HC0. After a maximum likelihood command such as poisson, expect Stata's robust standard error to be larger. The factor is the square root of n/(n - 1).

  • "The IRR from a modified Poisson model is an incidence rate ratio."

    Stata labels the column IRR because of the irr option. 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 and its time denominator, not from the option name.

What to do in your own analysis

Glossary

Poisson regression (การถดถอยแบบ Poisson)
A regression for counts that models the log of the mean count and assumes that the variance equals the mean.
empirical variance (ความแปรปรวนเชิงประจักษ์)
The average squared distance of the observations from their mean, dividing by n rather than n - 1.
overdispersion (การกระจายเกินกว่าที่แบบจำลองกำหนด)
More spread in a count outcome than the Poisson model assumes: the variance exceeds the mean.
underdispersion (การกระจายน้อยกว่าที่แบบจำลองกำหนด)
Less spread in a count outcome than the Poisson model assumes: the variance is below 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.
sandwich estimator (ตัวประมาณแบบแซนด์วิช)
A variance estimate that keeps the coefficient but reads the spread from the data: inverse bread, meat, inverse bread. Its standard error is the robust standard error.
HC0 (HC0)
The basic sandwich, with no small-sample factor; for a maximum likelihood command, Stata's vce(robust) multiplies its variance by n/(n - 1).
score contribution (score contribution ของแต่ละหน่วยสังเกต)
One observation's share of the estimating equation. In a Poisson model with a log link, it is the covariate values times the residual; with only an intercept, the residual itself.
bread (ส่วน bread ของ sandwich estimator)
The outer factor of the sandwich, J: how sharply the estimating equation changes with the coefficient; here n times the fitted mean.
meat (ส่วน meat ของ sandwich estimator)
The middle factor of the sandwich, K: the sum of the squared score contributions, read from the data.
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.
log-binomial model (แบบจำลอง log-binomial)
A regression for a yes-or-no outcome that models the log of the risk with the Bernoulli variance built in; it estimates a risk ratio directly.
modified Poisson regression (Poisson ดัดแปลง)
Poisson regression of a binary outcome with a sandwich standard error, used to estimate a risk ratio.
rate ratio (อัตราส่วนอัตราอุบัติการณ์)
The ratio of two mean counts per unit of time; Stata labels it IRR, short for incidence rate ratio.

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. StataCorp. Stata User's Guide: [U] 20.22. Stata documentation at stata.com (PDF). Accessed 2026-10-11. https://www.stata.com/manuals/u20.pdf
  4. 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
  5. 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

  • A robust standard error leaves the coefficient unchanged; only the standard error and the interval move.
  • For an intercept-only Poisson model, the bread is n times the fitted mean and the meat is the sum of squared residuals. The HC0 variance is the meat over the bread squared.
  • The sandwich can make the standard error smaller or larger: in the hand example it gives 0.071 and 0.400 against a model-based 0.224. Four patients show the arithmetic but are far too few to trust these standard errors.
  • For a maximum likelihood command, Stata's vce(robust) multiplies the HC0 variance by n/(n - 1), so label hand values HC0.
  • A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias.

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

0
Message for International and Thai ReadersUnderstanding My Medical Context in ThailandRead more →Message for International and Thai ReadersUnderstanding My Broader Content Beyond MedicineRead more →

Comments

No comments yet. Be the first to share your thoughts.

Sign in to comment