*! version 1.5.0 26aug2026

Title

svylet -- Wald omnibus F-test, Bonferroni pairwise comparisons,
and Compact Letter Display, for svy: mean/total/
proportion

Syntax

svylet
varname
[if] [in],
over(varname)
stat(statname)
[options]

where statname is one of mean, total, proportion, or
ratio.

Main
over(varname) variable that defines the groups to compare;
required

stat(statname) statistic to estimate and compare:
mean, total, proportion, or ratio; required

level(#) category of varname to test, when
stat(proportion); default is level(1). Ignored (with a note)
for mean, total, and ratio

denominator(varname) denominator variable, when
stat(ratio) -- varname (the main argument) is the numerator.
Required with stat(ratio), ignored (with a note) otherwise

alpha(#) significance level used for the Bonferroni pairwise
comparisons and the Compact Letter Display; default is alpha(0.05)

Vs-a-reference (optional)
ref(#) value of over() to use as a fixed baseline
category. Adds k-1 Bonferroni-corrected Wald contrasts (each other
category against ref()) alongside the all-pairs CLD -- a DIFFERENT
family of hypotheses, not a replacement; see
Remarks

Bootstrap (optional)
boot(#) number of bootstrap replications for a prepivoted
p-value of the omnibus F-test; default is boot(0), meaning the
analytic p-value only

bseed(#) random-number seed set with set seed before
the boot() replications; default lets Stata's current seed stand

varname must be numeric. fweights, pweights, and design
variables are not specified on svylet itself: the dataset must
already be svyset before calling svylet, exactly as for
svy:.

Description

svylet targets any complex-survey setting with a multistage,
stratified-cluster design and survey weights -- household budget and
labor-force surveys, demographic and health surveys, panel household
surveys, and similar public-use microdata -- where a researcher already
uses svy: and needs a design-based significance test across the
categories of a grouping variable (years, waves, regions, cohorts),
not just the point estimates that svy: itself reports.

svylet answers the question that svy: mean,
total, and proportion with over() leave open: are the
k group estimates actually different from one another, accounting for
the survey design? It runs the corresponding svy: command once
internally, then:

1. tests the omnibus null hypothesis that all k group
estimates are equal, using a Wald F-test with the Korn & Graubard (1990)
degrees-of-freedom adjustment (the same adjustment test/testparm
apply by default after svy:);

2. tests every pairwise comparison among the k groups, with a
Bonferroni correction for the k(k-1)/2 comparisons;

3. assigns a Compact Letter Display (CLD): groups that share at
least one letter are not significantly different from each other at
alpha(), after the Bonferroni correction.

Because svylet already ran svy: internally, it also exposes the
point estimates, standard errors, confidence limits, and sample sizes by
group through return, so a caller does not need to run the same
svy: command a second time just to build a point-estimate table next
to the test. See Stored results below, and see
tsvy for a command that uses exactly this to build a full
table (national/region/department, by year) in one pass.

Options

Main

over(varname) specifies the variable whose categories define the
groups to compare (for example, a year variable). It must have at least 2
distinct values in the estimation sample; svylet stops with an error
otherwise. varname must be NUMERIC -- svy: mean/total/
proportion/ratio do not accept a string over() at all
(svylet checks for this and stops with a clear message naming
encode as the fix, rather than letting the underlying svy:
call fail with its own less specific error). If varname has a value
label, svylet does not use it for the test itself (values are only
shown as returned by levelsof), so the mapping from a raw code to a
human-readable label, if you need one, is the caller's responsibility (see
the years() logic in tsvy for one way to do this without
depending on a value label).

stat(statname) chooses which svy: command svylet runs
internally: mean, total, proportion, or ratio
(added in version 1.2). Unlike svy: itself, svylet accepts
exactly one analysis variable per call (not a varlist, and not a
num/den expression for ratio -- see denominator()
below), because the test is defined over the categories of over()
for a single measured quantity. To test several variables, call
svylet once per variable.

level(#) selects, for stat(proportion) only, which category of
varname the test and point estimates refer to (the "success" value).
varname produces one row of output in svy: proportion for every
distinct value it takes; level(#) tells svylet which of those
rows to carry through the test. It has no effect for mean,
total, or ratio, and svylet prints a note if it is set
to something other than its default alongside any of those.

denominator(varname) names the denominator variable, required when
stat(ratio) -- the command's main argument, varname, is the
numerator. svylet builds the num/den expression that
svy: ratio expects internally; you never write the slash yourself.
e(b)/e(V) from svy: ratio ..., over() have the same
shape as mean/total (one equation, one column per
over() category), so the test, the Bonferroni comparisons, and the
Compact Letter Display work exactly the same way for ratio as for
mean/total. Ignored (with a note) for any other stat().

alpha(#) sets the significance level for the Bonferroni-adjusted
pairwise comparisons feeding the Compact Letter Display. It does not
affect the omnibus F-test's own p-value, only which groups end up sharing
a letter.

Bootstrap

boot(#) requests a bootstrap-calibrated p-value for the omnibus
F-test, in addition to the analytic one, using the prepivoting approach of
Beran (1988) as adapted for hypothesis testing by Hall & Wilson (1991):
pool the data from all over() groups, resample whole primary
sampling units (PSUs) with replacement within their original strata
(bsample, cluster() strata()), then randomly reassign
each resampled cluster as a whole to a pseudo-group, in the same
proportions (by cluster count) as the real over() groups, so that
the null hypothesis holds by construction in every replication. See
Remarks for why the reassignment happens by whole
cluster and not by individual observation, and References
for the supporting literature. boot(0), the default, skips this
entirely and reports only the analytic p-value.

boot() requires that the current svyset declare both a PSU
and a stratum (any of the internal names Stata's svy: may store
them under -- e(psu), e(psu1), e(su), e(su1), and
correspondingly for strata -- are tried); svylet stops with a clear
error, showing what each of those macros actually contained, if none of
them resolve to real variables.

bseed(#) passes a seed to set seed immediately before running
the boot() replications, so results are reproducible. If omitted,
whatever seed Stata's random-number generator is already on is used
as-is.

Remarks and examples

Remarks are presented under the following headings:

The omnibus F-test and the Korn-Graubard adjustment

Bonferroni pairs and the Compact Letter Display

ref(): comparing against a fixed baseline, not all pairs

Degenerate variance (proportions of exactly 0 or 1)

Why boot() resamples and reassigns by whole PSU

What svylet deliberately does not do

The omnibus F-test and the Korn-Graubard adjustment

The omnibus test is a Wald test of H0: all k group means/totals/
proportions are equal, built from the k-1 contrasts of every group
against the first one. This is mathematically identical regardless of
which group is used as the reference contrast -- any set of k-1
independent contrasts spanning the same space gives the same Wald
statistic, because it is a linear reparameterization of the same test.
(An earlier version of svylet accepted a refgroup() option
that appeared to let the caller change the reference group; it was removed
in version 1.1 once this invariance was confirmed algebraically -- there
was nothing for it to do. See the changelog at the top of svylet.ado.)

Stata's own test/testparm after a svy: command apply, by
default, the Korn & Graubard (1990) degrees-of-freedom adjustment: with
k-1 the dimension of the test and d = e(df_r) the design
degrees of freedom,

F = [(d - (k-1) + 1) / ((k-1)*d)] * W with
F ~ F(k-1, d-(k-1)+1)

where W is the raw Wald statistic. svylet applies the same
adjustment (confirmed by comparing svylet's output against
svy: regress + testparm on the same design and data).

Bonferroni pairs and the Compact Letter Display

Every one of the k(k-1)/2 pairs of over() categories is
compared with a design-based t-test, and the raw p-value is multiplied by
the number of pairs (capped at 1) -- the standard Bonferroni correction.
Groups are then assigned to the smallest number of letters such that two
groups share a letter if and only if they belong to some maximal subset of
groups that are all pairwise non-significant at alpha(). This is the
same logic used by pwcompare, cld and by the classic
letter-display output of ANOVA packages: shared letters mean "not
distinguishable here", not "equal".

ref(): comparing against a fixed baseline, not all pairs

GRUPO/CLD and ref() answer two DIFFERENT questions, and can
legitimately disagree on the same data -- neither is "wrong" when they do.
This matters in practice: comparing tsvy's CLD letters against a reference
script that only tested each year against the most recent one found
about 17% of the derived year-vs-2026 calls disagreeing, entirely because
the two procedures test different families of hypotheses (see the tsvy.ado
v1.5 changelog for the worked example).

GRUPO/CLD (the default, always computed) answers "which of these k
categories differ from EACH OTHER" -- all k(k-1)/2 pairs are
compared, and the Bonferroni correction multiplies each raw p-value by
k(k-1)/2. ref(#) answers a narrower, DIFFERENT question --
"which categories differ from THIS ONE baseline category" -- only k-1
comparisons are in that family, so the Bonferroni correction multiplies by
k-1 instead. A comparison that survives Bonferroni at k-1
comparisons can fail to survive it at k(k-1)/2 comparisons (and,
less intuitively but just as legitimately, the reverse can happen too, in
either direction), because the two corrections are controlling the
family-wise error rate over two DIFFERENT families of hypotheses, not the
same family with a different sample size. Dunn (1961) is the classical
reference for applying the Bonferroni inequality to multiple-comparison
problems, including comparing several groups against one control.

The k-1 comparisons ref() builds share the same baseline
category, so their test statistics are correlated with each other in a
known way -- an all-pairs procedure (this command's own CLD, Tukey's
method, or a plain Bonferroni split across k-1 comparisons as
ref() does) does not exploit that correlation and is therefore more
conservative than necessary for this specific family. Dunnett (1955, 1964)
derived a single-step procedure for exactly this "several treatments vs a
control" design that uses the correlation among the k-1 contrasts to
get sharper (less conservative, more powerful) critical values while still
controlling the family-wise error rate at the nominal level; see also Hsu
(1996, chapter 4) and Bretz, Hothorn & Westfall (2010, chapter 4) for
modern treatments and worked examples (the latter's R package
multcomp implements it as contrMat(..., type="Dunnett")).
ref() in this command uses the simpler Bonferroni split (Dunn 1961),
not the Dunnett single-step critical value -- it is always valid (Bonferroni
never inflates the false-positive rate, regardless of the correlation
structure) but somewhat conservative relative to Dunnett's for this
specific comparison family; implementing genuine Dunnett critical values
requires the multivariate t distribution, not attempted here. See
References below.

Degenerate variance (proportions of exactly 0 or 1)

A group with a proportion of exactly 0 or 1 -- common in small domains --
has an undefined or zero design-based variance. Rather than silently
treating that as "no difference" (a missing variance compared against a
number is treated as larger than any real number internally, which would
otherwise merge that group into whatever comparison happens to run last),
svylet flags those categories explicitly, marks their letter as
?, and excludes them from every candidate subset in the Compact
Letter Display. A warning naming the affected categories (by their
position within over(), not their raw value) is printed. Treat
? as "the test could not be computed here", never as "no
difference".

The omnibus F only drops the degenerate categories from its own
contrast -- it does not go missing just because one category out of
many is degenerate. As long as at least 2 categories have a defined,
positive variance, the omnibus F and its p-value are computed on that
subset (a note is printed naming how many of the over()
categories were used); with fewer than 2 usable categories there is
nothing left to test, so r(F_omnibus)/r(p_omnibus) are
returned missing. Point estimates and pairwise comparisons that do not
involve a degenerate category are unaffected either way.

Why boot() resamples and reassigns by whole PSU

boot()'s reassignment step must move whole primary sampling units
(clusters), never individual observations, from their original
over() group to a pseudo-group. Splitting a cluster's observations
across different pseudo-groups destroys the intra-cluster correlation that
the design-based variance estimator (and the observed F-statistic) assumes
when computing variance, understating the variability of the bootstrap
null-reference distribution and inflating the false-positive rate. A Monte
Carlo simulation under a true null hypothesis (10,000 simulated datasets,
199 bootstrap replications each) found the row-level reassignment used by
svylet 1.0 rejected a true null about 9.3% of the time at a nominal
5% level -- nearly double -- while the cluster-level reassignment used
from 1.1 onward rejected about 3.2% of the time (conservative, not
inflated). See AUDIT.md and sim/simulacion_bootstrap.py in the
svylet repository for the full simulation and its results, and
References below for the supporting literature
on why resampling and reassignment under a complex design must preserve
the cluster as the unit of exchangeability.

What svylet deliberately does not do

o svylet takes one analysis variable at a time, not a
varlist -- see tsvy for looping this over many variables
and many levels of aggregation at once.

o svylet does not cross over() with a second grouping
dimension in the same call (for example, "by year, separately within each
region"); subset the data with if and call svylet again for
each region, or use tsvy, which does exactly that
internally.

o the omnibus F-test result does not depend on, and cannot be
changed by, which over() category is treated as the reference (see
above).

Examples

The leading . before each line below is the command prompt, shown
here only because that is the standard convention for Stata help files
-- it is not part of the command. Copying a single-line example with the
. included works fine (Stata's do-file runner tolerates it), but
copying a multi-line foreach/forvalues block with a .
left on every line, including the body and the closing brace, breaks
Stata's parsing of the block. When copying a loop example into a do-file,
strip the leading . first.

Setup: a design with each observation as its own PSU (no real
clustering in auto.dta, this is only enough for the analytic path)

. sysuse auto, clear

. gen long psu_id = _n

. svyset psu_id

proportion: is the share of imported cars different across
repair-record categories?

. svylet foreign, over(rep78) stat(proportion) level(1)

mean: does average mpg differ across repair-record categories?

. svylet mpg, over(rep78) stat(mean)

total: does total vehicle weight differ across repair-record
categories?

. svylet weight, over(rep78) stat(total)

ratio: does the trunk-to-length ratio differ across
repair-record categories? denominator() is a separate option, not a
num/den expression typed into the main argument

. svylet trunk, over(rep78) stat(ratio) denominator(length)

over() must be numeric -- if the grouping variable you have is
a STRING (say, origin, holding "Domestic"/"Foreign" text), run it
through encode first, then pass the resulting numeric variable:

. decode foreign, generate(origin)

. encode origin, generate(origin_num)

. svylet mpg, over(origin_num) stat(mean)

With a bootstrap-calibrated p-value (needs a design with real PSUs
and strata; industry/south in nlsw88.dta are used here
only to have multi-row clusters to demonstrate the mechanism, not as a
real survey design):

. sysuse nlsw88, clear

. svyset industry, strata(south) singleunit(certainty)

. svylet wage, over(race) stat(mean) boot(200) bseed(20260824)

A looser significance level widens which groups end up sharing a
letter (fewer groups end up "significantly different"), without changing
the omnibus F-test itself, only the pairwise comparisons behind
GRUPO:

. svylet mpg, over(rep78) stat(mean) alpha(0.10)

ref(): compare every rep78 category against ONE fixed
baseline (here, 5 = "Excellent") instead of all pairs -- k-1
Bonferroni-adjusted contrasts (Dunn 1961), a DIFFERENT family of
hypotheses from the all-pairs CLD above (see
Remarks). This is the pattern for "did each
group change relative to a reference category/year", the question
tsvy's refyear() answers over time -- see its help for a
worked example with years:

. svylet mpg, over(rep78) stat(mean) ref(5)

. return list

. matrix list r(p_vsref)

Reading the per-category vs-baseline p-value programmatically
(missing at r(ref_idx)'s own position, and throughout if
ref() was not specified):

. forvalues i = 1/`r(k_categorias)' {

. if `i'' != `r(ref_idx)'' di "rep78 = " `r(nombre_categoria_`i'')' ///

. " vs baseline p = " el(r(p_vsref), `i'', 1)

. }

Reading the results back after the command runs -- the scalars and
matrices in Stored results above:

. svylet mpg, over(rep78) stat(mean)

. return list

. matrix list r(b)

Reading the per-category letters and estimates programmatically,
one category at a time (this is exactly the pattern tsvy
uses internally to build a table row by row):

. forvalues i = 1/`r(k_categorias)' {

. di "rep78 = " `r(nombre_categoria_`i'')' ///

. " estimate = " el(r(b), `i'', 1) ///

. " letter = " `r(letra_`i'')'

. }

Stored results

svylet stores the following in r():

Scalars

r(F_omnibus) omnibus Wald F-statistic (Korn-Graubard
adjusted); missing if not computable (see r(k_categorias) below)

r(p_omnibus) analytic p-value of the omnibus F-test

r(p_omnibus_boot) bootstrap-calibrated p-value; missing
unless boot() was specified and at least one replication succeeded

r(B_efectivo) number of boot() replications that
completed without error; missing if boot() was not specified

r(df_num) numerator degrees of freedom of the omnibus test
(k-1)

r(df_den) denominator degrees of freedom of the omnibus test
(Korn-Graubard adjusted)

r(df_raw) design degrees of freedom from e(df_r),
unadjusted -- use this one, not r(df_den), if building an
ordinary t-based confidence interval by hand

r(k_categorias) number of over() categories tested
(k)

r(ref_idx) position (1,...,k) of ref() within
over()'s categories; 0 if ref() was not specified

Macros

r(letra_i) Compact Letter Display code for the
ith category of over() (? if that category had
degenerate variance), i = 1,...,k

r(nombre_categoria_i) value of over() that
the ith category corresponds to, i = 1,...,k

Matrices

r(b) k x 1: point estimate for each over()
category, in the same order as r(nombre_categoria_i)

r(V) k x k: design-based variance-covariance
matrix of r(b)

r(n_ponderado) k x 1: weighted sample size
(e(_N_subp)) per category

r(n_sin_ponderar) k x 1: unweighted sample size
(e(_N)) per category

r(ci_lower) k x 1: lower confidence limit per category,
read directly from r(table) of the underlying svy: call (in
the logit-transformed scale Stata itself uses for proportion, not
reconstructed by hand)

r(ci_upper) k x 1: upper confidence limit per category,
same source as r(ci_lower)

r(p_vsref) k x 1: Bonferroni-adjusted (k-1
comparisons) p-value of each category against ref()'s category;
missing throughout if ref() was not specified, and at
r(ref_idx)'s own position always (a category is not tested against
itself)

r(p_vsref_raw) k x 1: same comparisons, raw (unadjusted)
p-value -- for audit/diagnostic use; apply your own correction if
r(p_vsref)'s plain Bonferroni is not what you want

r(b), r(V), r(n_ponderado), r(n_sin_ponderar),
r(ci_lower), and r(ci_upper) let a caller build a complete
point-estimate table (estimate, standard error via sqrt(diag(r(V))),
confidence interval, sample sizes) without running the underlying
svy: command a second time. tsvy is built entirely on
these returns.

References

Beran, R. 1988.
Prepivoting test statistics: A bootstrap view of asymptotic refinements.
Journal of the American Statistical Association 83(403): 687-697.

Bretz, F., T. Hothorn, and P. Westfall. 2010.
Multiple Comparisons Using R. Boca Raton, FL: CRC Press.

Canty, A. J., and A. C. Davison. 1999.
Resampling-based variance estimation for labour force surveys.
The Statistician 48(3): 379-391.

Davison, A. C., and D. V. Hinkley. 1997.
Bootstrap Methods and Their Application. Cambridge University Press.

Dunn, O. J. 1961.
Multiple comparisons among means.
Journal of the American Statistical Association 56(293): 52-64.

Dunnett, C. W. 1955.
A multiple comparison procedure for comparing several treatments with a
control.
Journal of the American Statistical Association 50(272): 1096-1121.

Dunnett, C. W. 1964.
New tables for multiple comparisons with a control.
Biometrics 20(3): 482-491.

Field, C. A., and A. H. Welsh. 2007.
Bootstrapping clustered data.
Journal of the Royal Statistical Society, Series B 69(3): 369-390.

Hall, P., and S. R. Wilson. 1991.
Two guidelines for bootstrap hypothesis testing.
Biometrics 47(2): 757-762.

Hsu, J. C. 1996.
Multiple Comparisons: Theory and Methods. Boca Raton, FL: Chapman &
Hall/CRC.

Korn, E. L., and B. I. Graubard. 1990.
Simultaneous testing of regression coefficients with complex survey data:
Use of Bonferroni t statistics.
The American Statistician 44(4): 270-276.

Rao, J. N. K., and C. F. J. Wu. 1988.
Resampling inference with complex survey data.
Journal of the American Statistical Association 83(401): 231-241.

Rust, K. F., and J. N. K. Rao. 1996.
Variance estimation for complex surveys using replication techniques.
Statistical Methods in Medical Research 5(3): 283-310.

Wolter, K. M. 2007.
Introduction to Variance Estimation. 2nd ed. Springer.

Author

Andres Talavera Cuya. Affiliation stated for identification purposes
only -- this software is not an official product of INEI and INEI bears
no responsibility for it. Developed independently by the author; views,
methods, and results are the author's own and do not necessarily reflect
the position of INEI. Distributed under the GNU General Public License v3
(https://www.gnu.org/licenses/gpl-3.0.txt).

Source, installation instructions, and the companion command
tsvy: https://github.com/atalaveracuya/svylet.
This is not (yet) an SSC package; download svylet.ado and this help
file into a directory on your adopath (or clone the repository
and add it with adopath ++ <path>).

Suggested citation: Talavera Cuya, A. 2026. svylet: Stata module to test
equality of survey-weighted means, totals, proportions, and ratios across
groups. Available from
https://github.com/atalaveracuya/svylet.

Also see

Online: tsvy, svy, svy postestimation,
mean, total, proportion, test,
testparm, pwcompare, bsample

En espanol: svylet_es