Which vce() When: Model-Based, Robust or Cluster-Robust Standard Errors
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
A regression coefficient and its standard error come from different calculations, and Stata's vce() option sets the second. In 600 simulated patients from 40 hospitals, half the hospitals run a discharge follow-up programme. The programme's incidence rate ratio for readmissions is 0.72 under every choice. The standard error of its logarithm is 0.0686 when the Poisson variance is trusted, 0.0816 with vce(robust) and 0.1721 with vce(cluster hospital). The cluster-robust version adds up each hospital's contributions to the estimate before squaring, so it allows correlation between patients of one hospital. Its 95% interval, 0.51 to 1.01, is the one that reflects a programme allocated by hospital. A negative binomial model gives almost the same ratio but still treats patients as independent. No variance estimator changes what is estimated or removes confounding by hospital. This article gives a decision table for vce(), the matching R calls, and Methods and Results sentences for count and binary outcomes.
Three standard errors for one programme
A hospital network's quality team is studying readmissions, the hospital admissions each patient has in the year after a discharge. Half of its hospitals run a discharge follow-up programme. The team wants to know whether their patients are readmitted less often. The data here are simulated: 600 patients, 15 from each of 40 hospitals.
The analysis plan says vce(robust). A reviewer asks for vce(cluster hospital), and a statistician suggests a negative binomial model instead, a count model that allows extra spread. The programme's incidence rate ratio, the ratio of mean admissions per patient in programme hospitals to that in the others, barely moves between the three. The confidence interval does: two versions exclude 1, and one includes it.
This article explains what vce() controls, what each choice assumes and how to choose from the design before any result is seen, then how to report it.
What vce() controls
A regression command estimates two things: the coefficients, and how much they would vary across repeated studies. The second is the variance-covariance matrix, which Stata calls the VCE. Each coefficient's estimated variance sits on its diagonal, and the covariance of each pair of coefficients sits off it. Standard errors are the square roots of the diagonal.
The vce() option chooses how this matrix is estimated [1]. It never touches the coefficients, so the point estimate stays where it was. It changes the standard errors, and with them every confidence interval and P value.
The covariances matter when a result combines coefficients. The variance of a sum of two coefficients, as lincom computes it, includes twice their covariance. nlcom uses the whole matrix for nonlinear combinations, so both inherit whatever vce() produced.
Stata: the model-based fit and its VCE
* ---- 1. Model-based standard errors: vce(oim), the default -------------------------
poisson admissions i.program i.older, irr vce(oim) nolog
estat vce. * ---- 1. Model-based standard errors: vce(oim), the default -------------------------
. poisson admissions i.program i.older, irr vce(oim) nolog
Poisson regression Number of obs = 600
LR chi2(2) = 68.92
Prob > chi2 = 0.0000
Log likelihood = -957.42519 Pseudo R2 = 0.0347
------------------------------------------------------------------------------
admissions | IRR Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7189098 .0493347 -4.81 0.000 .6284364 .8224084
1.older | 1.59639 .1112163 6.71 0.000 1.392637 1.829953
_cons | 1.301232 .080132 4.28 0.000 1.153285 1.468159
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.
. estat vce
Covariance matrix of coefficients of poisson model
| admissions
| 1. 1.
e(V) | program older _cons
-------------+-----------------------------------
admissions |
1.program | .00470929
1.older | -.00006111 .00485354
_cons | -.0019389 -.00296885 .00379229vce(oim), Stata's default, typed out. With irr, the table shows rate ratios and their standard errors on the ratio scale. The interval is computed on the log scale and then exponentiated, which is why it is not symmetric around the rate ratio. estat vce then prints the VCE of the coefficients on the log scale, where the diagonal entry for 1.program is the square of 0.0686, the log-scale standard error.R: the same fit and the same matrix
# ---- 1. Model-based standard errors (the inverse information, as Stata's vce(oim)) ----
fit <- glm(admissions ~ program + older, family = poisson, data = d)
print(summary(fit)$coefficients)
b <- coef(fit)
V_oim <- vcov(fit)
print(V_oim)> fit <- glm(admissions ~ program + older, family = poisson,
+ data = d)
> print(summary(fit)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.2633118 0.06158143 4.275832 1.904251e-05
program -0.3300193 0.06862416 -4.809084 1.516236e-06
older 0.4677449 0.06966719 6.713992 1.893708e-11
> b <- coef(fit)
> V_oim <- vcov(fit)
> print(V_oim)
(Intercept) program older
(Intercept) 0.003792272 -1.938898e-03 -2.968830e-03
program -0.001938898 4.709275e-03 -6.111248e-05
older -0.002968830 -6.111248e-05 4.853517e-03glm() is on the log scale, so the program row shows the standard error 0.0686 directly. vcov(fit) prints the same matrix as estat vce, with the intercept first.Model-based: trusting the model's variance
Poisson regression models a count: the logarithm of a patient's mean count is a linear function of the covariates, and the count's variance is assumed to equal its mean. Stata's default for it, vce(oim), inverts the observed information matrix, which measures how sharply the log-likelihood curves around the estimate. That curvature is computed as if the Poisson variance were right. The result is the model-based standard error: 0.0686 for the programme's coefficient, the log of its rate ratio.
Clinical counts often break the variance assumption. When counts vary more than their mean, called overdispersion, model-based standard errors tend to be too small. When they vary less, called underdispersion, they tend to be too large; Part 1 shows both directions with numbers.
A quick check is the Pearson dispersion statistic, which compares the squared residuals with the Poisson variance and is near 1 when that variance fits. Here it is 1.41, as the Stata output below shows. The extra spread has a known source in this simulation: patients of the same hospital share that hospital's influence on their admissions.
Stata: the Pearson dispersion statistic
* Pearson dispersion: sum of (y - mu)^2/mu over n - 3 (1 when the Poisson variance is right)
predict double mu, n
generate double pearson2 = (admissions - mu)^2/mu
quietly summarize pearson2
canon pearson_dispersion 4 r(sum)/(e(N) - 3). * Pearson dispersion: sum of (y - mu)^2/mu over n - 3 (1 when the Poisson variance is right)
. predict double mu, n
. generate double pearson2 = (admissions - mu)^2/mu
. quietly summarize pearson2
. canon pearson_dispersion 4 r(sum)/(e(N) - 3)
CANON g2.pearson_dispersion 1.4052predict stores each patient's fitted mean count. Each squared residual is divided by that mean, and the sum is divided by n - 3, the number of patients minus the three coefficients. quietly hides the summary table, so the statistic is printed only by the script's own CANON line. A small helper function in the script echoes each result; read the number at the end of the line, 1.41 when rounded.Robust: independent patients, variance in doubt
vce(robust) keeps the bread from the model but takes the meat from the data. Each patient's score contribution enters the meat as a square, $K = \sum_i s_i s_i^{\top}$, and the variance is the sandwich $J^{-1} K J^{-1}$. It stays valid in large samples when the variance assumption is wrong, provided the mean model is right and patients are independent [2].
Part 1 works the bread and the meat as numbers, including the factor $n/(n - 1)$ that Stata applies, where $n$ is the number of observations. The matrix form and its variants HC0 to HC3 are in the article on robust standard errors.
Here the robust standard error of the log rate ratio is 0.0816, 1.19 times the model-based value. The 95% interval for the rate ratio widens a little, to 0.61 to 0.84 from 0.63 to 0.82. It still treats the 600 patients as independent, which these data do not support.
Cluster-robust: independent hospitals, correlated patients
Patients of the same hospital share its staff, its discharge habits and its catchment area, so their counts tend to rise and fall together. The intraclass correlation (ICC) measures this: the correlation between the outcomes of two patients in the same cluster. The article on cluster randomized trials shows how it drives study design. A cluster-robust standard error allows any such correlation by changing only the meat.
$$S_g = \sum_{i \in g} s_i$$First, the score contributions of all patients in cluster $g$, here one hospital, are added into one cluster score, $S_g$.
$$K_{\text{cl}} = \sum_{g=1}^{G} S_g S_g^{\top}$$Then the meat is built from the $G$ cluster scores, one per hospital, where $G$ is the number of clusters. The variance becomes $J^{-1} K_{\text{cl}} J^{-1}$ [3]. Squaring a sum keeps the product of every pair of scores from the same hospital. When a hospital pushes all its patients' counts up, or all down, their residuals share a sign and these products enlarge the meat.
Two conditions come with it. The clusters must be independent of each other, and there must be enough of them, because the meat now rests on $G$ cluster scores rather than on $n$ patients [3]. Stata's poisson multiplies this variance by $G/(G - 1)$, and the R code asks vcovCL() for the same factor, so the two sets of results below agree.
Generalized estimating equations (GEE) use the same cluster-summed meat around a model that assumes a simple pattern of correlation between patients of one cluster (the working correlation, a guess at that pattern that need not be exactly right) [3]. The repeated-measures guide covers them.
Worked example: one programme, three standard errors
The simulated data were built so that the true values are known. There are 600 patients in 40 hospitals, 15 per hospital. The even-numbered hospitals run the programme, so 300 patients are in programme hospitals. A yes-or-no variable marks the 302 patients in an older age group.
Each patient's admissions over one year follow a Poisson distribution. Each hospital has its own factor, $e^{u}$, where $u$ is drawn from a normal distribution with mean 0 and standard deviation 0.4. The factor multiplies the mean of all that hospital's patients, so patients of one hospital resemble each other.
In a hospital without the programme and with a factor of 1, the mean is 1.2 for a younger patient. It is multiplied by 0.75 under the programme and by 1.6 for an older patient. Averaged over hospitals, every mean is a little higher, because the factor averages a little above 1.
The patients had 872 admissions in all, a mean of 1.45 per patient with a variance of 2.34. The mean was 1.22 in programme hospitals and 1.69 in the others. Every row of the table below fits the same model, Poisson regression of admissions on the programme and the age group; only vce() changes.
The programme under three variance estimators
| Option | What it assumes | SE, log rate ratio | 95% CI, rate ratio | SE relative to model-based |
|---|---|---|---|---|
vce(oim), model-based | Patients independent; variance equals the mean | 0.0686 | 0.63 to 0.82 | reference |
vce(robust) | Patients independent; any variance | 0.0816 | 0.61 to 0.84 | 1.19 |
vce(cluster hospital) | Hospitals independent; any correlation within a hospital | 0.1721 | 0.51 to 1.01 | 2.51 |
Reading the table
The rate ratio is 0.72 in every row, because vce() never moves the coefficient. The cluster-robust standard error is 2.11 times the robust one and 2.51 times the model-based one. Its interval, 0.51 to 1.01, includes 1, while the other two exclude it.
The reason is how the programme was allocated. Every patient in a hospital shares its programme status and its hospital factor, so the comparison rests on 40 hospitals, not on 600 independent patients. Only the cluster-robust standard error counts the information that way.
The age group behaves differently: each hospital has older and younger patients in similar proportions. Its shared factor raises or lowers both groups together, so it largely cancels from the age contrast. The age group's cluster-robust standard error, 0.0697, is no larger than the robust 0.0818; in this sample it is 0.85 times it. Clustering costs most for a covariate that is the same for everyone in a cluster.
Because the model has no hospital terms, 0.72 estimates a ratio averaged over hospitals. With a log link and a hospital factor unrelated to the programme, that averaged ratio equals the within-hospital ratio. The data were generated with that ratio set to 0.75. Only the constant absorbs the hospital factor, which is why the Stata panes show a constant above 1.2.
This sample is not a typical draw. It was picked from several simulated samples because its programme and comparison hospitals received similar hospital factors. Its rate ratio therefore sits closer to 0.75 than most samples of 40 hospitals would give, and all three intervals contain 0.75.
In other samples, the estimate moves by about as much as the cluster-robust standard error says. The two narrower intervals would then miss 0.75 far more often than their nominal 95% level allows. Coverage, how often such intervals contain the true value, is judged over many samples, as the article on robust standard errors does by simulation.
Stata: vce(robust)
* ---- 3. vce(robust): sandwich, each patient independent ----------------------------
poisson admissions i.program i.older, irr vce(robust) nolog. * ---- 3. vce(robust): sandwich, each patient independent ----------------------------
. poisson admissions i.program i.older, irr vce(robust) nolog
Poisson regression Number of obs = 600
Wald chi2(2) = 41.34
Prob > chi2 = 0.0000
Log pseudolikelihood = -957.42519 Pseudo R2 = 0.0347
------------------------------------------------------------------------------
| Robust
admissions | IRR std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7189098 .0586611 -4.04 0.000 .6126588 .8435876
1.older | 1.59639 .1306192 5.72 0.000 1.359854 1.87407
_cons | 1.301232 .0873629 3.92 0.000 1.140792 1.484237
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.Stata: vce(cluster hospital)
* ---- 4. vce(cluster hospital): sandwich, hospitals independent ----------------------
poisson admissions i.program i.older, irr vce(cluster hospital) nolog. * ---- 4. vce(cluster hospital): sandwich, hospitals independent ----------------------
. poisson admissions i.program i.older, irr vce(cluster hospital) nolog
Poisson regression Number of obs = 600
Wald chi2(2) = 45.15
Prob > chi2 = 0.0000
Log pseudolikelihood = -957.42519 Pseudo R2 = 0.0347
(Std. err. adjusted for 40 clusters in hospital)
------------------------------------------------------------------------------
| Robust
admissions | IRR std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7189098 .1237477 -1.92 0.055 .5130444 1.007381
1.older | 1.59639 .1112673 6.71 0.000 1.39255 1.830068
_cons | 1.301232 .156567 2.19 0.029 1.027866 1.647302
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.R: robust and cluster-robust standard errors
# ---- 3. Robust: sandwich, each patient independent (HC0 times n/(n - 1), as Stata's vce(robust)) ----
V_rob <- sandwich(fit) * n / (n - 1)
se_rob <- sqrt(diag(V_rob))
canon("robust.se_program", se_rob[["program"]])
irr_ci("robust.irr_program", b[["program"]], se_rob[["program"]])
canon("robust.se_older", se_rob[["older"]])
# ---- 4. Cluster-robust: hospitals independent (HC0 with the G/(G - 1) factor, as Stata's vce(cluster)) ----
V_cl <- vcovCL(fit, cluster = ~hospital, type = "HC0", cadjust = TRUE)
se_cl <- sqrt(diag(V_cl))
canon("cluster.n_clusters", length(unique(d$hospital)), 0)
canon("cluster.irr_program", exp(b[["program"]]))
canon("cluster.se_program", se_cl[["program"]])
irr_ci("cluster.irr_program", b[["program"]], se_cl[["program"]])
canon("cluster.se_older", se_cl[["older"]])
irr_ci("cluster.irr_older", b[["older"]], se_cl[["older"]])
canon("ratio.cluster_robust_program", se_cl[["program"]] / se_rob[["program"]])
canon("ratio.cluster_oim_program", se_cl[["program"]] / se_oim[["program"]])
canon("ratio.robust_oim_program", se_rob[["program"]] / se_oim[["program"]])
canon("ratio.cluster_robust_older", se_cl[["older"]] / se_rob[["older"]]) > V_rob <- sandwich(fit) * n/(n - 1)
> se_rob <- sqrt(diag(V_rob))
> canon("robust.se_program", se_rob[["program"]])
CANON g2.robust.se_program 0.0816
> irr_ci("robust.irr_program", b[["program"]], se_rob[["program"]])
CANON g2.robust.irr_program_lb 0.6127
CANON g2.robust.irr_program_ub 0.8436
> canon("robust.se_older", se_rob[["older"]])
CANON g2.robust.se_older 0.0818
> V_cl <- vcovCL(fit, cluster = ~hospital, type = "HC0",
+ cadjust = TRUE)
> se_cl <- sqrt(diag(V_cl))
> canon("cluster.n_clusters", length(unique(d$hospital)),
+ 0)
CANON g2.cluster.n_clusters 40
> canon("cluster.irr_program", exp(b[["program"]]))
CANON g2.cluster.irr_program 0.7189
> canon("cluster.se_program", se_cl[["program"]])
CANON g2.cluster.se_program 0.1721
> irr_ci("cluster.irr_program", b[["program"]], se_cl[["program"]])
CANON g2.cluster.irr_program_lb 0.5130
CANON g2.cluster.irr_program_ub 1.0074
> canon("cluster.se_older", se_cl[["older"]])
CANON g2.cluster.se_older 0.0697
> irr_ci("cluster.irr_older", b[["older"]], se_cl[["older"]])
CANON g2.cluster.irr_older_lb 1.3926
CANON g2.cluster.irr_older_ub 1.8301
> canon("ratio.cluster_robust_program", se_cl[["program"]]/se_rob[["program"]])
CANON g2.ratio.cluster_robust_program 2.1095
> canon("ratio.cluster_oim_program", se_cl[["program"]]/se_oim[["program"]])
CANON g2.ratio.cluster_oim_program 2.5083
> canon("ratio.robust_oim_program", se_rob[["program"]]/se_oim[["program"]])
CANON g2.ratio.robust_oim_program 1.1890
> canon("ratio.cluster_robust_older", se_cl[["older"]]/se_rob[["older"]])
CANON g2.ratio.cluster_robust_older 0.8518 sandwich(fit) is the uncorrected sandwich, multiplied by n/(n - 1) to match Stata's vce(robust). vcovCL() with cadjust = TRUE applies the G/(G - 1) factor of Stata's vce(cluster). Lines beginning CANON are printed by a small helper function in the script that echoes each result; read the number at the end of each line and ignore the prefix. They match the Stata output.Robust Poisson or negative binomial: two targets
Negative binomial regression keeps the Poisson mean model and adds a dispersion parameter, $\alpha$, so that a count with mean $\mu$ has variance $\mu + \alpha\mu^2$ [4]. It models the whole distribution of counts. Robust Poisson regression claims less: it models the mean and lets the sandwich handle the variance. The article on Poisson rates covers the model itself.
When the mean model is right, both estimate the same rate ratio, so the choice is about the target. Robust Poisson suits a question about the mean association alone. Negative binomial regression suits a question about the counts themselves, such as the share of patients with no admission. It can also be more precise when its variance model is close to the truth.
The two also weight patients differently. Negative binomial regression weights each patient's residual by $1/(1 + \alpha\mu)$, so a patient with a high expected count counts for less; Poisson regression weights all residuals alike.
If the mean model is wrong, for example if the programme works better in busier hospitals, the two can report different rate ratios. Each is then a different weighted average of the varying effects, and no choice of vce() reconciles them. The section What vce() cannot change, further down, comes back to this point.
In the simulated data, nbreg gives a rate ratio of 0.72, with a 95% interval of 0.62 to 0.84. It differs from the Poisson ratio only in the third decimal, because the weights differ. The estimated $\alpha$ is 0.23, and the likelihood-ratio statistic for $\alpha = 0$ is 34.41, so the counts are overdispersed. The standard error of the log rate ratio, 0.0795, is close to the robust Poisson value of 0.0816.
Both treat patients as independent. Here the overdispersion and the clustering share one source, the hospital factor, and the negative binomial model absorbs the extra variance without knowing who shares it. With vce(cluster hospital) added in Stata, its standard error rises to 0.1702 and its interval becomes 0.52 to 1.01. The last Stata output in this section shows that fit.
Stata: negative binomial regression
* ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
nbreg admissions i.program i.older, irr nolog. * ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
. nbreg admissions i.program i.older, irr nolog
Negative binomial regression Number of obs = 600
LR chi2(2) = 48.68
Dispersion: mean Prob > chi2 = 0.0000
Log likelihood = -940.21983 Pseudo R2 = 0.0252
------------------------------------------------------------------------------
admissions | IRR Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7226321 .0574734 -4.08 0.000 .6183271 .8445322
1.older | 1.590514 .1276263 5.78 0.000 1.35905 1.8614
_cons | 1.30083 .0918756 3.72 0.000 1.132665 1.493962
-------------+----------------------------------------------------------------
/lnalpha | -1.455862 .2286047 -1.903919 -1.007805
-------------+----------------------------------------------------------------
alpha | .2331992 .0533104 .1489836 .3650192
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to incidence-rate ratios.
Note: _cons estimates baseline incidence rate.
LR test of alpha=0: chibar2(01) = 34.41 Prob >= chibar2 = 0.000R: glm.nb from the MASS package
# ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
nb <- glm.nb(admissions ~ program + older, data = d)
print(summary(nb))> nb <- glm.nb(admissions ~ program + older, data = d)
> print(summary(nb))
Call:
glm.nb(formula = admissions ~ program + older, data = d, init.theta = 4.288179278,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.26300 0.07037 3.737 0.000186 ***
program -0.32486 0.07953 -4.085 4.41e-05 ***
older 0.46406 0.08024 5.783 7.32e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(4.2882) family taken to be 1)
Null deviance: 703.84 on 599 degrees of freedom
Residual deviance: 652.94 on 597 degrees of freedom
AIC: 1888.4
Number of Fisher Scoring iterations: 1
Theta: 4.288
Std. Err.: 0.980
2 x log-likelihood: -1880.440glm.nb() reports theta, which is 1/alpha, and log-scale coefficients that match the Stata rate ratios once exponentiated. The line saying the dispersion parameter is taken to be 1 refers to the GLM scale factor, not to alpha; the number inside Negative Binomial() is theta. The program standard error agrees with the Stata log-scale value to four decimals.Stata: negative binomial with vce(cluster hospital)
* negative binomial with a cluster-robust standard error (Stata only)
nbreg admissions i.program i.older, irr vce(cluster hospital) nolog
canon nb_cluster.se_program.stata 4 _se[1.program]. * negative binomial with a cluster-robust standard error (Stata only)
. nbreg admissions i.program i.older, irr vce(cluster hospital) nolog
Negative binomial regression Number of obs = 600
Wald chi2(2) = 44.64
Dispersion: mean Prob > chi2 = 0.0000
Log pseudolikelihood = -940.21983 Pseudo R2 = 0.0252
(Std. err. adjusted for 40 clusters in hospital)
------------------------------------------------------------------------------
| Robust
admissions | IRR std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7226321 .1229716 -1.91 0.056 .5176864 1.008713
1.older | 1.590514 .1109611 6.65 0.000 1.387248 1.823564
_cons | 1.30083 .1573644 2.17 0.030 1.026239 1.648892
-------------+----------------------------------------------------------------
/lnalpha | -1.455862 .3898306 -2.219916 -.6918084
-------------+----------------------------------------------------------------
alpha | .2331992 .0909082 .1086182 .5006699
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to incidence-rate ratios.
Note: _cons estimates baseline incidence rate.
. canon nb_cluster.se_program.stata 4 _se[1.program]
CANON g2.nb_cluster.se_program.stata 0.1702Choosing vce(): start from the design
The choice follows from how the data were collected, and it is best made before any standard error is seen. Two questions do most of the work: are the observations independent, and is the model's variance credible? If the observations are clustered, the clustering decides, whatever the answer to the second question. Cluster on the unit that was sampled or allocated, such as a hospital that ran a programme or a patient with repeated rows.
| Data structure | Stata | R | What to watch |
|---|---|---|---|
| Independent patients; variance credible, such as a dispersion statistic near 1 | vce(oim), the default | vcov(fit) | Rests entirely on the model's variance |
| Independent patients; counts over- or underdispersed | vce(robust) | sandwich(fit) * n / (n - 1) | The mean model must still be right |
| Independent patients; yes-or-no outcome, risk ratio wanted: modified Poisson regression, a Poisson model fitted to the 0 or 1 outcome | poisson y x, irr vce(robust) | glm(y ~ x, family = poisson) plus sandwich() | The sandwich is required, not optional |
| Patients clustered in hospitals, or repeated rows per patient | vce(cluster unit), such as vce(cluster hospital) or vce(cluster patient) | vcovCL(fit, cluster = ~unit, type = "HC0", cadjust = TRUE), such as ~hospital | Cluster on the unit sampled or allocated; clusters independent of each other. With patients nested in hospitals, cluster on the hospital, the outer unit. Or model the clustering (GEE, or a mixed model with its own term for each hospital), as the repeated-measures guide describes |
| Few clusters | vce(jackknife, cluster(hospital)), or a cluster-robust standard error judged against a t distribution | sandwich::vcovJK(fit, cluster = ~hospital) | The plain cluster-robust standard error tends to be too small |
| Resampling instead of a formula | vce(bootstrap, cluster(hospital)) or vce(jackknife, cluster(hospital)) | A bootstrap that resamples hospitals, or sandwich::vcovJK(fit, cluster = ~hospital) | Resample whole clusters, never single rows |
OPG, resampling and few clusters
One option looks like a sandwich but is not. vce(opg), the outer product of the gradients, inverts the meat alone, so like vce(oim) it is valid only when the model's variance is right. Here it gives 0.0595, smaller than even the model-based 0.0686. The Stata output after this section shows that fit.
The bootstrap and the jackknife refit the model on resampled data [1]. They respect clustering only when they resample whole hospitals, which in Stata happens only when cluster(hospital) is given inside vce() [1].
The cluster-robust standard error relies on many clusters. With few, it tends to be too small, and a small-sample correction, with significance judged against a t distribution, may then be preferable; the repeated-measures guide describes it. Report the number of clusters beside every cluster-robust result.
The 40 hospitals here, half of them with the programme, are usually treated as enough for the plain cluster-robust standard error. A small-sample correction would widen the interval somewhat and could be reported as a sensitivity analysis.
Stata: vce(opg)
* ---- 2. vce(opg): the outer product of the scores alone ---------------------------
poisson admissions i.program i.older, irr vce(opg) nolog
canon opg.se_program 4 _se[1.program]. * ---- 2. vce(opg): the outer product of the scores alone ---------------------------
. poisson admissions i.program i.older, irr vce(opg) nolog
Poisson regression Number of obs = 600
LR chi2(2) = 68.92
Prob > chi2 = 0.0000
Log likelihood = -957.42519 Pseudo R2 = 0.0347
------------------------------------------------------------------------------
| OPG
admissions | IRR std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.program | .7189098 .0427671 -5.55 0.000 .6397901 .807814
1.older | 1.59639 .0975423 7.66 0.000 1.416215 1.799488
_cons | 1.301232 .0741003 4.62 0.000 1.16381 1.454881
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.
. canon opg.se_program 4 _se[1.program]
CANON g2.opg.se_program 0.0595vce() option, the R call and the reporting sentences. It encodes the decision table above and computes no standard error.What vce() cannot change
Every option above changes how uncertain the estimate looks; none changes what is estimated. The estimand, the quantity the analysis sets out to estimate, is fixed by the question: the population, the outcome, the comparison and the summary measure. Here it is the ratio of mean admissions over one year in programme hospitals to that in other hospitals. It compares patients of the same age group, assumes one ratio for both age groups and averages over hospitals: a marginal, or population-averaged, rate ratio.
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.
Clustering does not remove confounding by hospital either. If hospitals that chose the programme also had healthier patients, the rate ratio would be biased under every vce(). In this simulation the programme followed the hospital number, independent of the hospital factor, so there is no such confounding.
Nor does a variance estimator change what a sample can estimate. Adding vce(robust) to a modified Poisson model fitted to a case-control sample does not turn its estimate into a risk ratio. How the controls were sampled decides the measure, as the article on control sampling shows.
Reporting: Methods and Results sentences
The Methods should name the variance estimator, the reason for it and the number of clusters. The Results should say which interval is reported.
For the yes-or-no outcome, any readmission within the year, the risk ratio comes from modified Poisson regression: a log-link Poisson model fitted to the 0 or 1 outcome, with a sandwich variance [5]. The sandwich is required because the Poisson variance, the risk $p$ itself, exceeds the Bernoulli variance of a yes-or-no outcome, $p(1 - p)$. Zou and Donner extended the method to clustered binary data with a cluster-robust variance [6]. The step-down that leads to this model is in the article on risk ratios for common outcomes.
In the simulated data, the risk of any readmission, printed in the first Stata output below, was 0.667 in programme hospitals and 0.763 in the others. Adjusted for age group, the risk ratio is 0.87, with a robust interval of 0.79 to 0.96. The cluster-robust interval, the one to report because the programme was allocated by hospital, is 0.75 to 1.01. The risk ratio is also averaged over hospitals and assumes one ratio for both age groups; the simulation set no true risk ratio.
Stata: the risk of any readmission in each group
quietly summarize any_admission if program == 1
canon any.risk_program 4 r(mean)
quietly summarize any_admission if program == 0
canon any.risk_control 4 r(mean). quietly summarize any_admission if program == 1
. canon any.risk_program 4 r(mean)
CANON g2.any.risk_program 0.6667
. quietly summarize any_admission if program == 0
. canon any.risk_control 4 r(mean)
CANON g2.any.risk_control 0.7633any_admission is 1 for a patient with at least one admission and 0 otherwise, so its mean in each group is the risk of any readmission. quietly hides the summary tables, so each risk is printed only by the script's own CANON line. A small helper function in the script echoes each result; read the number at the end of each line: 0.667 in programme hospitals and 0.763 in the others.Stata: modified Poisson, robust and cluster-robust
* ---- 6. Binary outcome (any admission): modified Poisson for the risk ratio ------------
poisson any_admission i.program i.older, irr vce(robust) nolog
canon any.rr_program 4 exp(_b[1.program])
canon any.se_robust 4 _se[1.program]
canon any.rr_lb_robust 4 exp(_b[1.program] - invnormal(0.975)*_se[1.program])
canon any.rr_ub_robust 4 exp(_b[1.program] + invnormal(0.975)*_se[1.program])
poisson any_admission i.program i.older, irr vce(cluster hospital) nolog. * ---- 6. Binary outcome (any admission): modified Poisson for the risk ratio ------------
. poisson any_admission i.program i.older, irr vce(robust) nolog
Poisson regression Number of obs = 600
Wald chi2(2) = 23.03
Prob > chi2 = 0.0000
Log pseudolikelihood = -569.62054 Pseudo R2 = 0.0058
-------------------------------------------------------------------------------
| Robust
any_admission | IRR std. err. z P>|z| [95% conf. interval]
--------------+----------------------------------------------------------------
1.program | .8709482 .0446817 -2.69 0.007 .7876326 .963077
1.older | 1.231833 .0645076 3.98 0.000 1.111672 1.364982
_cons | .6845157 .0322619 -8.04 0.000 .6241162 .7507604
-------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.
. canon any.rr_program 4 exp(_b[1.program])
CANON g2.any.rr_program 0.8709
. canon any.se_robust 4 _se[1.program]
CANON g2.any.se_robust 0.0513
. canon any.rr_lb_robust 4 exp(_b[1.program] - invnormal(0.975)*_se[1.program])
CANON g2.any.rr_lb_robust 0.7876
. canon any.rr_ub_robust 4 exp(_b[1.program] + invnormal(0.975)*_se[1.program])
CANON g2.any.rr_ub_robust 0.9631
. poisson any_admission i.program i.older, irr vce(cluster hospital) nolog
Poisson regression Number of obs = 600
Wald chi2(2) = 25.91
Prob > chi2 = 0.0000
Log pseudolikelihood = -569.62054 Pseudo R2 = 0.0058
(Std. err. adjusted for 40 clusters in hospital)
-------------------------------------------------------------------------------
| Robust
any_admission | IRR std. err. z P>|z| [95% conf. interval]
--------------+----------------------------------------------------------------
1.program | .8709482 .0657274 -1.83 0.067 .7511992 1.009787
1.older | 1.231833 .0558421 4.60 0.000 1.127106 1.346291
_cons | .6845157 .0400294 -6.48 0.000 .6103887 .7676449
-------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.vce(robust) and the second vce(cluster hospital); the ratios are identical and only the intervals change. Lines beginning CANON are printed by a small helper function in the script that echoes each result; read the number at the end of each line and ignore the prefix.Sentences to adapt
Each Methods sentence names the model, the measure, the variance estimator, the reason and the cluster unit. Each Results sentence names the interval it reports.
Methods, yes-or-no outcome: Risk ratios for any readmission within one year were estimated by modified Poisson regression adjusted for age group, with cluster-robust standard errors and the 40 hospitals as clusters, because the programme was allocated by hospital.
Results, yes-or-no outcome: Any readmission occurred in 66.7% of patients in programme hospitals and 76.3% in other hospitals (adjusted risk ratio 0.87, 95% CI 0.75 to 1.01, from cluster-robust standard errors).
Methods, count outcome: Incidence rate ratios for admissions in the year after discharge were estimated by Poisson regression adjusted for age group, with cluster-robust standard errors and the 40 hospitals as clusters, because the programme was allocated by hospital. Negative binomial regression with the same clustering was a sensitivity analysis.
Results, count outcome: Patients in programme hospitals had a mean of 1.22 admissions against 1.69 in other hospitals (adjusted incidence rate ratio 0.72, 95% CI 0.51 to 1.01, from cluster-robust standard errors; negative binomial with the same clustering 0.72, 95% CI 0.52 to 1.01).
Follow-up here is one year for every patient. If it varied, the count model would need log follow-up time as an offset, giving a rate per unit of time, as the article on Poisson rates explains.
Common misreadings and their fixes
-
"vce(cluster) fixes confounding by hospital."
Here the rate ratio is 0.72 under every
vce()option, so any bias by hospital would survive all of them.Fix: Clustering changes the standard error only. Confounding by hospital needs the model, for example hospital as a covariate or a design that removes it. For an exposure set per hospital, as here, hospital indicators would absorb it, so measured hospital traits or randomising hospitals must do that work.
-
"Robust standard errors fix a misspecified model."
The sandwich guards the standard error against a wrong variance assumption, nothing more. A mean model that misses a curve or an interaction answers a different question.
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.
-
"Pick the vce() whose P value looks best."
Here the model-based interval excludes 1 and the cluster-robust interval includes it, so choosing after looking would choose the conclusion.
Fix: Write the variance estimator into the analysis plan from the design: the cluster unit, the reason and the expected number of clusters.
-
"A negative binomial model takes care of clustering."
It models extra variance but still treats patients as independent: its standard error here is 0.0795, against 0.1721 for cluster-robust Poisson.
Fix: Choose negative binomial or robust Poisson for the target, then cluster the standard error on the hospital in either model.
-
"Any identifier will do as the cluster variable."
Clustering on a patient identifier when every patient has one row returns the robust standard error and ignores the hospitals. Too few clusters tend to give a standard error that is too small.
Fix: Cluster on the unit that was sampled or allocated, report the number of clusters, and use a small-sample correction when they are few.
What to do in your own analysis
- Name the unit that was sampled or allocated, and write the
vce()choice and its reason into the analysis plan before fitting. - Check the dispersion of counts, but let the design, not the check, decide whether to cluster.
- Report the number of clusters with every cluster-robust result, and consider a small-sample correction when they are few.
- State whether the target is the mean association or the count distribution before choosing robust Poisson or negative binomial regression.
- Run
lincomornlcomafter the fit for any contrast that combines coefficients, so that the chosen VCE carries through.
Glossary
- variance-covariance matrix (VCE) (เมทริกซ์ความแปรปรวนและความแปรปรวนร่วม)
- 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.
- Poisson regression
- A model for counts in which the log of the mean count is a linear function of the covariates and the variance is assumed to equal the mean.
- incidence rate ratio
- The ratio of two event rates; with the same follow-up for everyone, the ratio of mean counts per person.
- observed information matrix (OIM)
- The curvature of the log-likelihood at the estimate; its inverse is the model-based VCE that vce(oim) reports.
- 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 formula that places the meat between two inverted copies of the bread, so the variance comes from the data rather than from the model's variance assumption.
- robust (sandwich) standard error (ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช)
- The standard error from the sandwich estimator; in large samples it stays valid when the variance assumption is wrong, if the mean model is right and observations are independent.
- bread (ส่วน bread ของ sandwich estimator)
- The outer factor of the sandwich, written J: how sharply the estimating equations change with the coefficients; it appears inverted on both sides of the meat.
- meat (ส่วน meat ของ sandwich estimator)
- The middle factor of the sandwich, written K: the summed squares of the score contributions, read from the data.
- score contribution (score contribution ของแต่ละหน่วยสังเกต)
- One observation's share of the estimating equation, written s_i; in Poisson regression, its covariate values times its residual.
- 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.
- cluster-robust standard error
- A sandwich standard error whose meat adds up the score contributions within each cluster before squaring, so observations in one cluster may be correlated; it needs many independent clusters.
- intraclass correlation (ICC)
- The correlation between the outcomes of two members of the same cluster.
- generalized estimating equations (GEE) (สมการประมาณค่าวางนัยทั่วไป)
- A method for correlated data that fits a mean model that assumes a simple pattern of correlation between members of one cluster (the working correlation, a guess at that pattern that need not be exactly right) and uses the cluster-summed sandwich for its standard errors.
- coverage
- How often intervals built the same way contain the true value over many samples; a 95% interval should do so 95% of the time.
- negative binomial regression
- A count model with the Poisson mean structure plus a dispersion parameter, so the variance can exceed the mean.
- outer product of the gradients (OPG)
- A variance estimate that inverts the meat alone; like the model-based one, it is valid only when the model's variance is right.
- estimand
- The quantity an analysis sets out to estimate, fixed by the question: the population, the outcome, the comparison and the summary measure. The variance estimator does not change it.
- marginal (population-averaged) rate ratio
- A rate ratio averaged over the clusters, here the hospitals, rather than one that compares patients within the same hospital.
- modified Poisson regression (Poisson ดัดแปลง)
- Poisson regression with a log link fitted to a yes-or-no outcome, with a sandwich standard error, to estimate a risk ratio.
- Bernoulli variance (ความแปรปรวนแบบ Bernoulli)
- The variance of a yes-or-no outcome with probability p, p(1 - p); it is below the Poisson value p whenever p lies between 0 and 1.
References
- StataCorp. Stata Base Reference Manual: [R] vce_option. Stata documentation at stata.com (PDF); release and year not stated. Accessed 2026-10-11. https://www.stata.com/manuals/rvce_option.pdf
- StataCorp. Stata User's Guide: [U] 20.22, Obtaining robust variance estimates. Stata documentation at stata.com (PDF); release and year not stated. Accessed 2026-10-11. https://www.stata.com/manuals/u20.pdf
- Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73:13-22. doi:10.1093/biomet/73.1.13 https://doi.org/10.1093/biomet/73.1.13
- Gardner W, Mulvey EP, Shaw EC. Regression analyses of counts and rates: Poisson, overdispersed Poisson, and negative binomial models. Psychol Bull. 1995;118(3):392-404. doi:10.1037/0033-2909.118.3.392 https://doi.org/10.1037/0033-2909.118.3.392
- 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
- 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
Key takeaways
- The vce() option changes the variance-covariance matrix and every standard error, interval and P value built on it, never the coefficient: the programme's rate ratio is 0.72 under every choice.
- Cluster on the unit that was sampled or allocated; with the programme allocated by hospital, the cluster-robust standard error was 2.51 times the model-based one, and its interval included 1.
- Robust Poisson and negative binomial regression share a mean model but differ in target and in how they weight patients, and neither handles clustering until its standard error is clustered.
- 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.
- Choose the variance estimator from the design before seeing results, and report the model, the measure, the estimator and the number of clusters in one Methods sentence.
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]]