Singleton PSUs in Survey Analysis with Python

Detecting and resolving strata with a single primary sampling unit

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 and compare strategies: certainty, skip, center, scale, collapse, and pool.

Author

Mamadou S. Diallo, Ph.D.

Published

July 21, 2026

Modified

July 21, 2026

Keywords

singleton PSU Python, lonely PSU survey, single PSU stratum, survey variance estimation Python, stratum collapse Python, lonely.psu 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 resolve them before you can estimate standard errors. This tutorial shows how svy detects singletons through the sample.singleton interface and walks through every handling strategy.

A Sample with Singletons

We’ll work with a small stratified design so the structure stays visible. It has five regions; East and West were each covered by a single cluster, making them singletons.

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,
    "cluster": ["c01", "c01", "c02", "c02", "c03", "c03", "c04", "c04",
                "c05", "c05", "c06", "c06", "c07", "c07", "c07", "c08", "c08"],
    "hhsize":  [4, 2, 3, 5, 6, 3, 2, 4, 5, 3, 4, 2, 3, 6, 2, 4, 1],
    "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"),
)

The sample.singleton accessor is the single entry point for everything that follows — detection, diagnostics, and handling.

Detecting Singletons

The quick checks answer “do I have a problem, and how big is it?”

print("Any singletons? ", sample.singleton.exists)
print("How many?       ", sample.singleton.count)
print("Which strata?   ", sample.singleton.keys())

Any singletons?  True
How many?        2
Which strata?    ['East', 'West']

show() returns a tidy table of the offending strata — the singleton key, how many records it holds, and its lone PSU:

print(sample.singleton.show())

shape: (2, 4)
┌───────────────┬───────┬─────┬────────┐
│ singleton_key ┆ n_obs ┆ psu ┆ region │
│ ---           ┆ ---   ┆ --- ┆ ---    │
│ str           ┆ i64   ┆ str ┆ str    │
╞═══════════════╪═══════╪═════╪════════╡
│ East          ┆ 3     ┆ c07 ┆ East   │
│ West          ┆ 2     ┆ c08 ┆ West   │
└───────────────┴───────┴─────┴────────┘

summary() adds context and, importantly, a recommended strategy based on how prevalent the singletons are:

info = sample.singleton.summary()

print(f"Singletons: {info.n_singletons} of {info.n_strata} strata")
print(f"Rows affected: {info.pct_rows_affected:.1f}%")
print(f"Recommendation: {info.recommendation}")
print(f"Why: {info.recommendation_reason}")

Singletons: 2 of 5 strata
Rows affected: 29.4%
Recommendation: SingletonHandling.POOL
Why: Multiple singletons (2); pooling creates a valid pseudo-stratum for variance estimation

To see every stratum side by side — PSU counts, sizes, and which ones are singletons — use strata_profile(). Pass a variable to include its (weighted) mean per stratum:

print(sample.singleton.strata_profile("income"))

shape: (5, 5)
┌─────────────┬────────┬───────┬─────────────┬──────────────┐
│ stratum_key ┆ n_psus ┆ n_obs ┆ mean_income ┆ is_singleton │
│ ---         ┆ ---    ┆ ---   ┆ ---         ┆ ---          │
│ cat         ┆ u32    ┆ u32   ┆ f64         ┆ bool         │
╞═════════════╪════════╪═══════╪═════════════╪══════════════╡
│ Center      ┆ 2      ┆ 4     ┆ 1147.5      ┆ false        │
│ East        ┆ 1      ┆ 3     ┆ 906.666667  ┆ true         │
│ North       ┆ 2      ┆ 4     ┆ 1030.0      ┆ false        │
│ South       ┆ 2      ┆ 4     ┆ 1060.0      ┆ false        │
│ West        ┆ 1      ┆ 2     ┆ 790.0       ┆ true         │
└─────────────┴────────┴───────┴─────────────┴──────────────┘

Why Singletons Block Variance Estimation

By default, svy refuses to estimate a variance while unhandled singletons remain, 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.singleton.summary(), then pick a strategy:    │
│                                                                         │
│    • sample.singleton.certainty()  — treat as self-representing units   │
│    • sample.singleton.skip()       — drop from variance (R 'remove')    │
│    • sample.singleton.center()     — grand-mean centering (R 'adjust')  │
│    • sample.singleton.scale()      — variance inflation (R 'average')   │
│    • sample.singleton.collapse()   — merge into nearby strata           │
│    • sample.singleton.pool()       — pool singletons into one stratum   │
│    • sample.singleton.combine(map) — manual stratum/PSU remapping       │
│                                                                         │
│  where estimation                                                       │
│                                                                         │
╰─────────────────────────────────────────────────────────────────────────╯

You resolve this by choosing a handling strategy. The rest of the tutorial covers each one.

Handling Strategies

svy offers a full menu, spanning the same ground as R’s lonely.psu option plus a few more:

Strategy	Method	Idea	R `lonely.psu`
Certainty	`certainty()`	Treat each singleton PSU as self-representing (no first-stage variance)	—
Skip	`skip()`	Exclude singleton strata from the variance calculation	`"remove"`
Center	`center()`	Score the singleton by its squared deviation from the grand mean	`"adjust"`
Scale	`scale()`	Drop singletons, then inflate variance by 1 / (1 − singleton fraction)	`"average"`
Collapse	`collapse()`	Merge each singleton into a nearby non-singleton stratum	—
Pool	`pool()`	Combine all singletons into a single pseudo-stratum	—
Combine	`combine(mapping)`	Manual stratum/PSU remapping	—

Every method returns a new Sample; the original is untouched. Applying one makes estimation work:

handled = sample.singleton.pool()
print(handled.estimation.mean("income"))

╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮
│                                                         │
│         est        se        lci          uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  1,014.7059   37.6010   910.3087   1,119.1030     3.71  │
│                                                         │
╰─────────────────────────────────────────────────────────╯

Handling records a recipe — it doesn’t rewrite your strata

svy never modifies your original design columns. Instead, a handling method attaches a variance recipe that the estimator reads at estimation time, and records what it did in last_result. That means re-inspecting the handled sample still reports the original singletons (handled.singleton.count is unchanged) — that is expected. To see the effect, look at last_result:

r = handled.singleton.last_result
print(f"{r.method}: {r.n_strata_before} strata → {r.n_strata_after}")

SingletonHandling.POOL: 5 strata → 4

Collapse: merge into nearby strata

collapse() is the most common choice when singletons are few. It merges each singleton into an existing non-singleton stratum. Before committing, inspect the merge targets for a singleton with candidates_for():

print(sample.singleton.candidates_for("East"))

shape: (3, 3)
┌─────────────┬────────┬───────┐
│ stratum_key ┆ n_psus ┆ n_obs │
│ ---         ┆ ---    ┆ ---   │
│ str         ┆ i64    ┆ i64   │
╞═════════════╪════════╪═══════╡
│ Center      ┆ 2      ┆ 4     │
│ North       ┆ 2      ┆ 4     │
│ South       ┆ 2      ┆ 4     │
└─────────────┴────────┴───────┘

By default candidates are ranked by size (smallest first), so using="smallest" merges into the smallest available stratum. You can also let svy match on similarity of one or more variables:

mapping = sample.singleton.suggest_mapping(variables=["income", "hhsize"])
print("Suggested merges:", mapping)

Suggested merges: {'East': 'North', 'West': 'North'}

Apply the collapse. Use within= to constrain merges to strata sharing a value (e.g. only merge inside the same super-region), and using= to pick the strategy — "smallest", "largest", "next"/"previous" by an ordering, an explicit {singleton: target} dict, or a callable:

collapsed = sample.singleton.collapse(using="smallest")

r = collapsed.singleton.last_result
print(f"Collapsed: {r.n_strata_before} strata → {r.n_strata_after}")
print(collapsed.estimation.mean("income"))

Collapsed: 5 strata → 3

╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮
│                                                         │
│         est        se        lci          uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  1,014.7059   55.9667   870.8388   1,158.5730     5.52  │
│                                                         │
╰─────────────────────────────────────────────────────────╯

Pool: gather singletons into one pseudo-stratum

When singletons don’t naturally belong to any existing stratum, pool() combines them all into a single new pseudo-stratum with enough PSUs for variance estimation, leaving the well-formed strata alone:

pooled = sample.singleton.pool(name="__pooled__")
print(pooled.estimation.mean("income"))

╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮
│                                                         │
│         est        se        lci          uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  1,014.7059   37.6010   910.3087   1,119.1030     3.71  │
│                                                         │
╰─────────────────────────────────────────────────────────╯

Certainty, skip, scale, and center

These four keep the sample structure and differ only in how the singleton contributes to the variance. certainty() treats each singleton PSU as self-representing; skip() drops it from the variance (R’s "remove"); scale() drops it but inflates the remaining variance (R’s "average"); and center() scores it by its deviation from the grand mean (R’s "adjust"):

for name, smp in {
    "certainty": sample.singleton.certainty(),
    "skip":      sample.singleton.skip(),
    "scale":     sample.singleton.scale(),
    "center":    sample.singleton.center(),
}.items():
    print(f"{name:10s} applied →", smp.singleton.last_result.method)

certainty  applied → SingletonHandling.CERTAINTY
skip       applied → SingletonHandling.SKIP
scale      applied → SingletonHandling.SCALE
center     applied → SingletonHandling.CENTER

Comparing the strategies

Because each method returns a Sample, you can run the same estimate through all of them and compare the standard errors. The point estimate barely moves; the strategies differ in how conservative the variance is:

strategies = {
    "certainty": sample.singleton.certainty(),
    "skip":      sample.singleton.skip(),
    "scale":     sample.singleton.scale(),
    "center":    sample.singleton.center(),
    "collapse":  sample.singleton.collapse(using="smallest"),
    "pool":      sample.singleton.pool(),
}

comparison = pl.concat(
    smp.estimation.mean("income").to_polars().select(
        pl.lit(name).alias("strategy"),
        pl.col("est").round(1),
        pl.col("se").round(2),
        pl.col("cv").round(4),
    )
    for name, smp in strategies.items()
)

print(comparison)

shape: (6, 4)
┌───────────┬────────┬───────┬────────┐
│ strategy  ┆ est    ┆ se    ┆ cv     │
│ ---       ┆ ---    ┆ ---   ┆ ---    │
│ str       ┆ f64    ┆ f64   ┆ f64    │
╞═══════════╪════════╪═══════╪════════╡
│ certainty ┆ 1014.7 ┆ 40.05 ┆ 0.0395 │
│ skip      ┆ 1079.2 ┆ 52.23 ┆ 0.0484 │
│ scale     ┆ 1014.7 ┆ 40.46 ┆ 0.0399 │
│ center    ┆ 1014.7 ┆ 49.21 ┆ 0.0485 │
│ collapse  ┆ 1014.7 ┆ 55.97 ┆ 0.0552 │
│ pool      ┆ 1014.7 ┆ 37.6  ┆ 0.0371 │
└───────────┴────────┴───────┴────────┘

The `handle()` dispatcher

If the strategy is chosen at runtime (e.g. read from a config), handle() dispatches by name and forwards any extra arguments:

strategy = "collapse"  # could come from a config file

result = sample.singleton.handle(strategy, using="smallest")
print(result.estimation.mean("income"))

╭──────────────── Estimate: MEAN (TAYLOR) ────────────────╮
│                                                         │
│         est        se        lci          uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  1,014.7059   55.9667   870.8388   1,158.5730     5.52  │
│                                                         │
╰─────────────────────────────────────────────────────────╯

Choosing a Strategy

There is no universally correct choice — it depends on how many singletons you have and why:

Few singletons, similar strata nearby — collapse() (optionally within= a super-region) preserves the design most faithfully.
Several singletons, no natural home — pool() builds one valid pseudo-stratum.
Singletons affect very few rows — skip() has minimal impact.
You want a conservative, automatic adjustment — scale() or center() mirror the familiar R behaviors.
The singletons are genuine certainty units — certainty() is the principled choice.

When in doubt, start from sample.singleton.summary().recommendation, which encodes these heuristics.

Next Steps

With singletons resolved, your design produces valid standard errors and you can move on to more specialized analyses.

Ready for more analysis?
Continue to Categorical Data Analysis →