import polars as pl
import svy
from rich import print as rprint
survey = pl.DataFrame({
"region": ["North"]*4 + ["South"]*4 + ["Center"]*4 + ["East"]*3 + ["West"]*2,
"zone": ["N"]*4 + ["S"]*4 + ["N"]*4 + ["N"]*3 + ["S"]*2,
"cluster": ["c01", "c01", "c02", "c02", "c03", "c03", "c04", "c04",
"c05", "c05", "c06", "c06", "c07", "c07", "c07", "c08", "c08"],
"income": [980, 1220, 1050, 870, 760, 1180, 990, 1310,
1420, 1090, 880, 1200, 950, 1010, 760, 720, 860],
"w": [20.0]*17,
})
sample = svy.Sample(
survey,
svy.Design(stratum="region", psu="cluster", wgt="w"),
)Singleton PSUs in Survey Analysis with Python
Detecting strata with a single primary sampling unit and declaring how to handle them
Tutorials
Survey Estimation
Variance Estimation
Python
Detect and handle singleton (lonely) PSUs in complex survey designs with svy. Learn why a stratum with a single PSU breaks variance estimation, declare a rule on the design with svy.Singleton, and compare strategies: center, scale, skip, self-representing, collapse, and pool.
Keywords
singleton PSU Python, lonely PSU survey, single PSU stratum, survey variance estimation Python, stratum collapse Python, lonely.psu Python, survey.adjust.domain.lonely Python, singleunit Stata Python, complex survey design Python, stratified sampling variance Python
A singleton PSU (also called a lonely PSU) is a stratum that contains only one primary sampling unit. Singletons are a common obstacle in complex survey analysis: variance estimation for a stratified multistage design measures the variability between PSUs within each stratum, and that requires at least two PSUs per stratum. A stratum with a single PSU offers no within-stratum variability to measure, so the variance is undefined.
Singletons arise for several reasons: fine stratification, certainty units, nonresponse emptying a stratum, or restricting a full sample to a domain. Whatever the cause, you have to say how to handle them before you can estimate standard errors. In svy you do that once, on the design, with svy.Singleton, as you would with R’s lonely.psu option or Stata’s singleunit(). The sample reports which strata are singletons, like it reports its strata and PSUs.
A Sample with Singletons
We’ll work with a small stratified design so the structure stays visible. It has five regions in two zones; East and West were each covered by a single cluster, making them singletons.
Detecting Singletons
sample.singletons lists the strata with a single PSU, one row each, with the lone PSU and the stratum’s rows; sample.n_singletons counts them. They sit next to the sample’s other design facts, sample.n_strata and sample.n_psus:
print(f"{sample.n_singletons} singletons among {sample.n_strata} strata and {sample.n_psus} PSUs")
print(sample.singletons)2 singletons among 5 strata and 8 PSUs
shape: (2, 4)
┌────────┬─────────┬─────┬─────────┐
│ region ┆ cluster ┆ n ┆ handled │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ u32 ┆ str │
╞════════╪═════════╪═════╪═════════╡
│ East ┆ c07 ┆ 3 ┆ null │
│ West ┆ c08 ┆ 2 ┆ null │
└────────┴─────────┴─────┴─────────┘The handled column is empty for now: the design has no singleton rule yet.
Why Singletons Block Variance Estimation
By default, svy refuses to estimate a variance while singletons remain without a rule, rather than silently returning a wrong standard error. Attempting an estimate raises a SingletonError that names the offending strata and lists your options:
from svy.errors import SvyError
try:
sample.estimation.mean("income")
except SvyError as err:
rprint(err)╭──────────────────────────────────────────────────────────────────────────────────────╮ │ │ │ ✗ 2 singleton PSU(s) detected [SINGLETON_ERROR] │ │ │ │ Found 2 singleton PSU(s) in the following strata: │ │ 1. region=East (PSU=c07, n=3) │ │ 2. region=West (PSU=c08, n=2) │ │ │ │ Variance cannot be estimated with unhandled singleton PSUs. │ │ Inspect them with sample.singletons, then declare a rule: │ │ │ │ • svy.Singleton("center") — grand-mean centering (R "adjust") │ │ • svy.Singleton("scale") — variance scaled up (R "average") │ │ • svy.Singleton("skip") — contribute nothing (R "remove") │ │ • svy.Singleton("collapse") — merge into another stratum │ │ • svy.Singleton("pool") — pool the singletons into one stratum │ │ • svy.Singleton("self_representing") — the PSU becomes a stratum │ │ │ │ where estimation │ │ hint sample.update_design(singleton=svy.Singleton("center")), or svy.Design(..., │ │ singleton="center") when building the sample. │ │ │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
Declaring a Rule
A singleton rule is part of the design. Pass it to svy.Design, as a method name or as svy.Singleton(...) when it takes options:
centered = svy.Sample(
survey,
svy.Design(stratum="region", psu="cluster", wgt="w", singleton="center"),
)
print(centered.estimation.mean("income"))╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮ │ y: income │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 1,014.7059 49.2126 858.0894 1,171.3223 4.85 │ │ │ ╰─────────────────────────────────────────────────────────╯
On an existing sample, update_design declares or replaces it, and records the change in design_history like any design edit:
pooled = sample.clone().update_design(singleton=svy.Singleton("pool"))
print(pooled.design.singleton)Singleton('pool')method |
Idea | R lonely.psu |
Stata singleunit() |
|---|---|---|---|
"center" |
Score the singleton by its squared deviation from the grand mean | "adjust" |
centered |
"scale" |
Leave singleton strata out, then inflate the variance by 1 / (1 − singleton fraction) | "average" |
scaled |
"skip" |
Singleton strata contribute nothing to the variance | "remove" |
certainty |
"self_representing" |
The singleton’s PSU becomes a stratum and its SSUs (or rows) the PSUs | — | — |
"collapse" |
Merge each singleton into another stratum | — | — |
"pool" |
Combine all singletons into one pseudo-stratum | — | — |
Note that R’s and Stata’s “certainty” means a singleton contributes nothing to the variance: that is svy’s "skip". "self_representing" is different: it keeps the singleton’s variance, measured one stage down.
NoteA rule is intent: it follows your data
svy never modifies your design columns. The rule tells the estimator what to do with whatever strata have a single PSU, and the handled column of sample.singletons shows what it did on the current data:
print(pooled.singletons)shape: (2, 4)
┌────────┬─────────┬─────┬────────────────────┐
│ region ┆ cluster ┆ n ┆ handled │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ u32 ┆ str │
╞════════╪═════════╪═════╪════════════════════╡
│ East ┆ c07 ┆ 3 ┆ pool -> __pooled__ │
│ West ┆ c08 ┆ 2 ┆ pool -> __pooled__ │
└────────┴─────────┴─────┴────────────────────┘The rule is applied again whenever the data or the design changes. Drop cluster c02, and North becomes a singleton too; the declared rule handles it, with nothing to redo:
fewer = pooled.wrangling.filter_records(svy.col("cluster") != "c02")
print(fewer.singletons)
print(fewer.estimation.mean("income"))shape: (3, 4)
┌────────┬─────────┬─────┬────────────────────┐
│ region ┆ cluster ┆ n ┆ handled │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ u32 ┆ str │
╞════════╪═════════╪═════╪════════════════════╡
│ East ┆ c07 ┆ 3 ┆ pool -> __pooled__ │
│ North ┆ c01 ┆ 2 ┆ pool -> __pooled__ │
│ West ┆ c08 ┆ 2 ┆ pool -> __pooled__ │
└────────┴─────────┴─────┴────────────────────┘╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮ │ y: income │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 1,022.0000 53.3193 873.9618 1,170.0382 5.22 │ │ │ ╰─────────────────────────────────────────────────────────╯
A sample with no singletons keeps the rule, idle, until it has some. The rule is saved and restored with the design, and to_code() writes it:
print(pooled.to_code())import svy
sample = svy.Sample(
data,
svy.Design(
stratum='region',
wgt='w',
psu='cluster',
singleton=svy.Singleton('pool'),
),
)
Collapse: merge into nearby strata
collapse is the most common choice when singletons are few. It merges each singleton into an existing stratum with several PSUs. By default the target is the smallest available stratum (using="smallest"). within= restricts the targets to strata with the same value of a column constant within each stratum, here the zone:
collapsed = svy.Sample(
survey,
svy.Design(
stratum="region", psu="cluster", wgt="w",
singleton=svy.Singleton("collapse", within="zone"),
),
)
print(collapsed.singletons)
print(collapsed.estimation.mean("income"))shape: (2, 4)
┌────────┬─────────┬─────┬────────────────────┐
│ region ┆ cluster ┆ n ┆ handled │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ u32 ┆ str │
╞════════╪═════════╪═════╪════════════════════╡
│ East ┆ c07 ┆ 3 ┆ collapse -> Center │
│ West ┆ c08 ┆ 2 ┆ collapse -> South │
└────────┴─────────┴─────┴────────────────────┘╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮ │ y: income │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 1,014.7059 57.4435 867.0428 1,162.3690 5.66 │ │ │ ╰─────────────────────────────────────────────────────────╯
using= also takes "largest", "next"/"previous" (by order_by=), a callable, or an explicit mapping {singleton: target} in the stratum column’s values:
explicit = sample.clone().update_design(
singleton=svy.Singleton("collapse", using={"East": "North", "West": "South"})
)
print(explicit.singletons["handled"].to_list())['collapse -> North', 'collapse -> South']An explicit mapping names strata by their values, so it can stop fitting the data after a filter or recode. svy then raises at the next estimate, naming the strata to add to the mapping; entries for strata that are no longer singletons are simply ignored.
Pool: gather singletons into one pseudo-stratum
When singletons don’t naturally belong to any existing stratum, pool combines them all into a new pseudo-stratum with enough PSUs for variance estimation, leaving the well-formed strata alone:
pooled = sample.clone().update_design(singleton=svy.Singleton("pool", name="others"))
print(pooled.estimation.mean("income"))╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮ │ y: income │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 1,014.7059 37.6010 910.3087 1,119.1030 3.71 │ │ │ ╰─────────────────────────────────────────────────────────╯
Comparing the methods
The same estimate under every method shows how much the choice matters. The point estimate does not move; the methods differ in how conservative the variance is:
rules = {
"center": "center",
"scale": "scale",
"skip": "skip",
"self_representing": "self_representing",
"collapse": svy.Singleton("collapse", within="zone"),
"pool": "pool",
}
comparison = pl.concat(
sample.clone()
.update_design(singleton=rule)
.estimation.mean("income")
.to_polars()
.select(
pl.lit(name).alias("method"),
pl.col("est").round(1),
pl.col("se").round(2),
pl.col("cv").round(4),
)
for name, rule in rules.items()
)
print(comparison)shape: (6, 4)
┌───────────────────┬────────┬───────┬────────┐
│ method ┆ est ┆ se ┆ cv │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ f64 ┆ f64 ┆ f64 │
╞═══════════════════╪════════╪═══════╪════════╡
│ center ┆ 1014.7 ┆ 49.21 ┆ 0.0485 │
│ scale ┆ 1014.7 ┆ 47.6 ┆ 0.0469 │
│ skip ┆ 1014.7 ┆ 36.87 ┆ 0.0363 │
│ self_representing ┆ 1014.7 ┆ 40.05 ┆ 0.0395 │
│ collapse ┆ 1014.7 ┆ 57.44 ┆ 0.0566 │
│ pool ┆ 1014.7 ┆ 37.6 ┆ 0.0371 │
└───────────────────┴────────┴───────┴────────┘Because the rule is a value, a method chosen at runtime (read from a config, say) is just svy.Singleton(method, **options).
Singletons Inside a Domain
A stratum can hold several PSUs in the sample but only one with rows in an estimation domain: a where= subpopulation or a by= level. Take households with an income under 1,000. Center has two clusters, but only c06 has such a household. sample.domain_singletons() lists these strata:
low_income = svy.col("income") < 1000
print(centered.domain_singletons(where=low_income))shape: (1, 4)
┌────────┬─────────┬─────┬────────┐
│ region ┆ cluster ┆ n ┆ n_psus │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ str ┆ u32 ┆ u32 │
╞════════╪═════════╪═════╪════════╡
│ Center ┆ c06 ┆ 1 ┆ 2 │
└────────┴─────────┴─────┴────────┘Two settings of the rule decide what happens. domains= picks the variance:
"standard"(default) uses the usual domain variance: the PSUs without domain rows count as zero totals, so the variance is well defined. It is R’s default."apply"(centerandscaleonly) handles such a stratum as a singleton, R’soptions(survey.adjust.domain.lonely = TRUE). Undercenterall the stratum’s PSU totals (the domain’s and the zero totals of the others) are centered at the grand mean; underscalethe stratum is left out and counted with the singletons.
on_domain_singletons= says what svy tells you:
on_domain_singletons= |
Reported |
|---|---|
"ignore" (default) |
The finding is kept on the result and in sample.warnings; nothing is printed |
"warn" |
Also a note under the result |
"error" |
Raises before estimating, naming the strata |
The default prints nothing because in designs with two PSUs per stratum (NHANES, MEPS, NHIS) almost every domain table has such strata. The finding is still there for anyone who checks:
result = centered.estimation.mean("income", where=low_income)
print(result.findings[0].code, result.findings[0].extra["pairs"])DOMAIN_SINGLETON_PSU ['region=Center']With on_domain_singletons="warn" the note appears under the table, and with domains="apply" the note says what was adjusted:
applied = sample.clone().update_design(
singleton=svy.Singleton("center", domains="apply", on_domain_singletons="warn")
)
print(applied.estimation.mean("income", where=low_income))╭─────────────────────────────────────── Estimate: MEAN (TAYLOR) ────────────────────────────────────────╮ │ y: income │ │ where: income < 1000 │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 863.3333 32.8379 722.0433 1,004.6233 3.80 │ │ │ │ note: 1 stratum with one PSU in the domain centered at the grand mean (domains="apply"): region=Center │ ╰────────────────────────────────────────────────────────────────────────────────────────────────────────╯
on_domain_singletons="error" stops instead, which is useful in a pipeline that must not publish these estimates unreviewed:
strict = sample.clone().update_design(
singleton=svy.Singleton("center", on_domain_singletons="error")
)
try:
strict.estimation.mean("income", where=low_income)
except SvyError as err:
rprint(err)╭──────────────────────────────────────────────────────────────────────────────────────╮ │ │ │ ✗ Strata with a single PSU in a domain [DOMAIN_SINGLETON] │ │ │ │ 1 stratum has one PSU in the domain: │ │ 1. region=Center │ │ │ │ The singleton rule was declared with domains="standard" and │ │ on_domain_singletons="error". │ │ Declare how to report them instead: │ │ │ │ • svy.Singleton("center", domains="apply") — treat them as singletons (R │ │ survey.adjust.domain.lonely = TRUE) │ │ • svy.Singleton("center", on_domain_singletons="warn") — note under the result │ │ • svy.Singleton("center") — recorded only (R default) │ │ │ │ where estimation │ │ param on_domain_singletons │ │ hint sample.domain_singletons(by=..., where=...) lists them. │ │ │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
NoteWhich rows count as the domain’s
A PSU is in the domain when it has a row inside where= and the by= level with the analysis variables present, whatever that row’s weight. This matches R, whose subset() keeps zero-weight rows. The degrees of freedom and n still count only rows with a nonzero weight, as R’s do. A calibrated design (after poststratify, rake, calibrate) never has such strata: its scores are nonzero outside the domain, and R keeps every row.
Replicate Weights
Replicate weights built from a design with singletons follow the rule. Under collapse, pool and self_representing, create_jk_wgts, create_bs_wgts and create_brr_wgts use the rule’s strata and PSUs, so the replicates and the Taylor variance agree. Under skip, the singletons contribute nothing, which is what the jackknife does with them. center and scale adjust the Taylor variance and have no replicate analogue; the jackknife and the bootstrap then raise, suggesting collapse:
jk = collapsed.weighting.create_jk_wgts()
print(jk.estimation.mean("income", method="replication"))╭────────────── Estimate: MEAN (JACKKNIFE) ───────────────╮ │ y: income │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 1,014.7059 57.9021 865.8638 1,163.5480 5.71 │ │ │ ╰─────────────────────────────────────────────────────────╯
Choosing a Method
There is no universally correct choice; it depends on how many singletons you have and why:
- Matching R or Stata defaults for published estimates:
center(R"adjust", Statacentered), the usual choice for DHS, NHANES and MEPS. - Few singletons, similar strata nearby:
collapse, optionallywithin=a region, preserves the design most faithfully. - Several singletons, no natural home:
poolbuilds one valid pseudo-stratum. - Singletons affect very few rows:
skiphas minimal impact. - The singletons’ PSUs have subsampled units:
self_representingmeasures their variance one stage down. - Domain estimates that must match R with
survey.adjust.domain.lonely = TRUE:svy.Singleton("center", domains="apply")or"scale".
TipComing from sample.singleton.center()?
Earlier versions handled singletons with methods on a sample.singleton accessor. It is gone: the rule is part of the design, so there is one place to set it. sample.singleton.center() is now svy.Design(..., singleton="center") or sample.update_design(singleton="center"); .certainty() is "self_representing"; .collapse(using=...) and .combine(mapping) on strata are svy.Singleton("collapse", using={singleton: target}), and .combine() on PSUs is sample.wrangling.recode(psu_column, {new: [old, ...]}, replace=True). The inspection methods became sample.singletons, sample.n_singletons and sample.domain_singletons(). Accessing sample.singleton raises an error with these spellings.
Next Steps
With singletons handled, your design produces valid standard errors and you can move on to more specialized analyses.
Ready for more analysis?
Continue to Categorical Data Analysis →