Understanding the VCE: Robust Standard Errors, the Sandwich Estimator and vce(cluster id)

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

Abstract

A regression coefficient and its uncertainty come from two separate calculations, and Stata's vce() option changes only the second. A standard error is not the spread of the data: with variance 18 across 100 hospitals, the standard error of the mean is 0.4243. Stata estimates a variance-covariance matrix, the VCE; standard errors are the square roots of its diagonal, and its covariances are needed for any sum of coefficients. The model-based VCE trusts the model's variance, while the robust sandwich puts the meat, built from squared residuals, between two inverted copies of the bread. Worked by hand for 100 hospitals, an intercept-only Poisson model gives a model-based standard error of 0.0447 and a robust one of 0.0849, with the same coefficient. In modified Poisson regression for a binary outcome, the robust standard error can be the smaller one. This article then contrasts vce(robust) with vce(cluster id), lists what a robust standard error cannot repair, and gives Methods and Results sentences.

🔍 Click to enlarge Click anywhere to close

Same coefficient, different uncertainty

A research team counts the patients admitted per day in 100 hospitals and fits a Poisson model in Stata with vce(robust), because the last analysis did. A second dataset follows the same patients over several admissions, and its code says vce(cluster patient_id). Many analysts read both options as commands that simply adjust the standard error, without knowing why the adjustment is needed or what happens inside the model.

This article builds the answer from the ground up: variance and the standard error, the variance-covariance estimator (VCE), the model-based and the robust standard error, and the sandwich estimator, then their use in Poisson regression, modified Poisson regression and clustered data.

Why a standard error exists

Suppose we study the number of patients admitted per day in 100 hospitals. Once the data are in, the mean is 5 admissions per day:

A single mean says nothing about how much the hospitals differ from each other.

Take two datasets with the same mean. In dataset A the counts spread little: mean 5, variance 1. In dataset B they spread a lot: mean 5, variance 16. Figure 1, set out as numbers in the table below, shows eight hospitals from each so that the pattern is visible.

If we collected data from another 100 hospitals, the mean estimated from a population like B would change more from sample to sample than one from a population like A, assuming the units are independent and identically distributed within each population. That is the reason a standard error exists.

Figure 1: same mean, different spread

Variance and standard error are different things

A variance describes how spread out the data are. A standard error (SE) describes how uncertain an estimate is over repeated sampling. With $n$ independent observations that each have variance $\sigma^2$, the variance of the mean is

so its standard error is

Suppose the admission counts of 100 hospitals have mean 5 and variance 18. Then

The estimated mean carries a sampling uncertainty of about 0.4243, on the same scale as the admission counts. The value 0.4243 is not the standard deviation of the hospitals; it is the standard error of the estimated mean.

Now keep the variance at 18 but raise the number of hospitals from 100 to 400:

The standard error falls by about half, because the number of observations rose fourfold. The sample size and the spread of the data both drive the standard error, and the two have to be kept apart.

What a VCE is

VCE stands for variance-covariance estimator: the method used to estimate the variances and covariances of the statistical estimates. When we write in Stata

poisson admissions treatment age, vce(robust)

the option vce(robust) means: use the robust variance-covariance estimator to estimate the uncertainty of the regression coefficients [1].

What Stata computes is not one standard error. It is the variance-covariance matrix of all the regression coefficients.

Why it has to be a matrix

Take a regression model

where $\beta_0$ is the intercept and $\beta_1$ the regression coefficient of $X$. The analysis gives the estimates

Each estimate has its own uncertainty, and the two may also change together when the data are sampled again. The VCE therefore has to estimate a matrix of this shape:

Where the VCE lives in Stata

After a regression model is fitted, Stata stores the main results in the objects listed in the table below. They can be displayed with

matrix list e(b)
matrix list e(V)
estat vce

If the model has 4 coefficients including the intercept, e(b) is 1 × 4 and e(V) is 4 × 4. This is the standard structure of Stata estimation results [1].

How the VCE changes the confidence interval and the P value

Suppose a regression gives $\hat\beta_1 = 0.40$ with a standard error of 0.10. A Wald 95% confidence interval on the coefficient scale, by the normal approximation, is

Now switch to a robust VCE. The standard error becomes 0.19, while the coefficient is still 0.40:

Figure 2: the same coefficient, two intervals

The model-based VCE

Any regression software needs a way to estimate the variance of the coefficients. The standard method in maximum likelihood estimation is the model-based VCE, which uses the likelihood and the model's own assumptions; its square roots are the model-based standard errors. In many models it is estimated as the inverse of the observed information matrix, the negative Hessian:

where

This is the formula of the observed-information VCE, not a formula that can stand in for every kind of VCE.

To make it intuitive, think of the Hessian as telling how quickly the likelihood changes when the regression coefficient changes. If the likelihood is strongly curved around the estimate, the coefficient is usually pinned down fairly precisely. If the likelihood is rather flat, the estimate of the coefficient is more uncertain.

The example of Poisson regression

Suppose we fit a Poisson regression model,

with

or, equivalently, $\mu_i = \exp(X_i^{\top}\boldsymbol\beta)$. The Poisson distribution has one key property:

The conditional variance equals the conditional mean. If the model expects a hospital to have 5 admissions on average, $E(Y_i \mid X_i) = 5$, the Poisson model also uses $Var(Y_i \mid X_i) = 5$ when it computes the model-based uncertainty.

What if the real variance is 18?

Now suppose the conditional variance is higher than the mean the Poisson assumption implies:

This is called overdispersion. When the Poisson model uses a variance of 5 but the real variability is higher, the model-based standard error can understate the uncertainty.

One caution. Comparing the mean and the variance of the whole dataset does not yet prove conditional overdispersion when the model has several predictors, because differences in the mean between groups also raise the overall variance. What has to be assessed is the variance conditional on the predictors in the model.

The robust VCE

A robust VCE estimates the variance-covariance matrix while relying less on some of the model's assumptions, above all the assumption about the variance. The most widely used version is the Huber-White sandwich variance estimator, and its square roots are the robust (sandwich) standard errors [2].

In Poisson regression, if the conditional mean model is right but the conditional variance does not follow the Poisson distribution, the regression coefficients can still be estimated with the Poisson estimating equations, and the robust VCE can be used for statistical inference. This idea is the basis of Poisson pseudo-maximum likelihood and of modified Poisson regression.

Model-based and robust VCE compared

How the sandwich estimator is computed

Once it is clear that the variance of the coefficients has to be estimated, the next step is how. The sandwich estimator has the basic formula

This is the basic form for independent observations, without any finite-sample correction. It is called a sandwich because it has three parts:

• Bread, $A^{-1}$, on the left
• Meat, $B$, in the middle
• Bread, $A^{-1}$, on the right

What the bread is

For Poisson regression with a log link,

where $X_i$ is the vector of predictors, including the intercept. The bread links the variability of the estimating equations to the variability of the estimated coefficients: it reflects how sensitive the equations used to estimate a coefficient are to a change in that coefficient.

What the meat is

First comes the score contribution of each observation. For Poisson regression it is

where $Y_i - \hat\mu_i$ is the raw residual, the difference between the observed value and the value the model predicts. The score contributions of all observations then build the meat:

Substituting the Poisson score gives

The formula uses each observation's squared residual. When the data vary a lot, the score contributions carry that variability into the meat.

Why the model-based and robust VCE differ

If the Poisson variance assumption is right,

so in expectation the meat agrees with the bread:

The sandwich VCE is then close to

which is the model-based VCE. But if the conditional variance does not follow the Poisson assumption,

the meat and the bread no longer need to be equal, and so

This is the mathematical reason the robust standard error differs from the model-based one. Stata uses the Huber-White sandwich in vce(robust) and extends the same idea to the cluster-robust VCE.

Worked example: the robust standard error by hand

Reproducing the hand example in Stata

The same dataset can be built and analysed in Stata without any real patient data. These are the hand example's commands:

* Create synthetic admission count data
clear

input admissions frequency
0  20
1  25
5  10
9  25
10 20
end

* Expand to 100 independent observations
expand frequency
drop frequency

* Examine the data
count
summarize admissions, detail

* Model 1: Conventional Poisson
poisson admissions
estimates store poisson_model

display "Coefficient = " _b[_cons]
display "Model-based SE = " _se[_cons]

* Model 2: Robust Poisson
poisson admissions, vce(robust)
estimates store poisson_robust

display "Coefficient = " _b[_cons]
display "Robust SE = " _se[_cons]

* Inspect coefficients and VCE
matrix list e(b)
matrix list e(V)
estat vce

* Compare the two models
estimates table poisson_model poisson_robust, b se

After the run, look for the results in the table below.

Stata: the 100 hospitals, model-based and robust

    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

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
------------------------------------------------------------------------------

Modified Poisson regression: why it needs the robust standard error

A second application matters a great deal in clinical epidemiology: estimating a risk ratio (RR) from a binary outcome with modified Poisson regression, a log-link Poisson model fitted to a 0 or 1 outcome. Suppose we study the association between inadequate antenatal care (ANC) and neonatal morbidity within 28 days. The outcome is binary:

What variance does a binary outcome have?

For a binary outcome with event probability $p$,

but its variance, the Bernoulli variance, is

unlike the Poisson model, which sets $Var(Y) = p$. Take a risk of 0.40:

A risk ratio from binary data

The Stata command for modified Poisson

* Binary outcome: 0 = No morbidity, 1 = Morbidity
* Independent observations

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

Stata heads the table IRR, its label for an incidence rate ratio, because the command is poisson, irr. But when the outcome is binary and the model estimates the risk within one follow-up horizon, with no person-time offset, $\exp(\hat\beta)$ is read as a risk ratio under a log-risk mean model. The IRR label Stata prints does not turn the estimate into an incidence rate ratio.

Adding exposure(tmorbidity28) to use person-time as the denominator changes the mean model into a rate model, and the effect measure has to be interpreted anew.

One more caution: a log-link risk model can predict probabilities above 1 for some covariate profiles. A robust standard error does not fix that.

vce(robust) or vce(cluster patient_id)?

Up to here every observation was assumed independent of the others. Many kinds of clinical data are not. Suppose we study the hospital admissions of 100 patients, each followed over several periods. Observations from the same patient may be related to each other.

Plain vce(robust) treats every row as an independent observation, so it may misjudge the uncertainty. The cluster-robust VCE, and the cluster-robust standard error it gives, is what has to be considered.

How the cluster-robust VCE works

In the ordinary robust VCE the meat uses the score contribution of each observation:

When there are several patients and each has several rows, the score contributions within the same patient are added together first. Let $g$ be a cluster, here a patient:

The cluster-level meat is then built from these sums:

which gives

with the finite-sample corrections not shown. This lets observations within a cluster be correlated, under two key conditions: the clusters must be independent of each other, and there must in general be enough independent clusters [5]. The same approach has been applied to modified Poisson regression for correlated binary outcomes [4].

The Stata code

* Independent observations
poisson admissions i.treatment c.age, vce(robust)

* Repeated observations within each patient
poisson admissions i.treatment c.age, ///
    vce(cluster patient_id)

vce(cluster patient_id) does not make the software build matched pairs, fix confounding or add patient fixed effects. It changes how the variance is estimated, by using the cluster as the independent unit when the score variability is computed.

With few clusters, the standard cluster-robust standard error may give confidence intervals and P values that cannot be trusted. Small-sample adjustments, or a method of inference suited to the number of clusters, then have to be considered.

The kinds of VCE in Stata

Stata supports several kinds of VCE; which ones are available depends on the estimation command [1].

OPG, bootstrap and jackknife

vce(opg) is not the robust sandwich VCE, even though it also uses score information, because the robust sandwich uses both the bread and the meat.

vce(bootstrap) and vce(jackknife) are not better than a robust standard error in every case. With clustered data, the resampling unit has to match the design of the study.

Decision table for clinical research

What a robust standard error cannot do

This point needs to be clear, because a robust standard error is easy to use for the wrong purpose. The table sets out which problems it can and cannot solve.

A robust standard error does not make the model right

Suppose the true conditional mean is

but we fit a model without $X^2$:

Even with

poisson outcome x, vce(robust)

we still cannot conclude that the coefficient of $X$ is the right value for the relationship in the true model. The robust VCE may give the uncertainty of a pseudo-true coefficient, the value a misspecified model converges to, but that does not automatically make the coefficient the clinical parameter or the causal effect we want.

In the same way, fitting modified Poisson regression to a case-control sample, in which cases and controls were sampled in proportions that do not reflect the source population, and reporting an RR without regard to the sampling design is not fixed by adding vce(robust). The article on control sampling shows how the design decides the measure.

Robust Poisson or negative binomial?

For a count outcome with overdispersion, both robust Poisson and negative binomial regression can be used, but the two do not have the same aim in every case.

Robust Poisson sets the conditional mean

and uses the sandwich VCE, without needing the Poisson variance to be right. If the conditional mean is right and the independence assumptions are suitable, Poisson pseudo-maximum likelihood gives consistent coefficient estimates even when the conditional variance differs from the conditional mean.

Negative binomial regression is another model, one that sets an explicit mean-variance relationship, such as the NB2 form

where $\alpha$ is the overdispersion parameter [6]. If the negative binomial distribution and variance structure fit, the model may be more efficient for some parameters, and it is useful for simulating or predicting the count distribution. That does not make it better than robust Poisson every time.

Choose from the aim of the study, the adequacy of the mean model and the distribution assumptions, not from the finding of overdispersion alone. The article on Poisson rates covers the count models themselves.

Reporting the results

Methods, modified Poisson regression:

Adjusted risk ratios and 95% confidence intervals were estimated using modified Poisson regression with a log link and robust sandwich variance estimation. The analysis assumed independent observations and a correctly specified conditional mean model.

If the data are clustered within hospitals and a cluster-robust VCE was actually used, the sentence changes to:

Adjusted risk ratios and 95% confidence intervals were estimated using modified Poisson regression with a log link and cluster-robust sandwich variance estimation, accounting for within-hospital correlation.

The two sentences mean different things and should not be swapped without checking the study design.

Results: suppose the analysis gave $aRR = 2.00$, 95% CI 1.26 to 3.17. It can be reported as:

Inadequate antenatal care was associated with a higher risk of neonatal morbidity within 28 days (adjusted RR 2.00, 95% CI 1.26 to 3.17).

These numbers are an example for writing only, not the result of analysing a real dataset. A manuscript must use the numbers of the adjusted model that was actually run.

Summary: what has to be understood

Related in the wiki: [[sandwich-standard-error-worked-by-hand]] [[heteroskedasticity-robust-standard-errors]] [[repeated-measures-modeling-guide]] [[poisson-regression-rates-lambda]] [[cluster-randomized-trials-design-ethics]]