Cluster Sampling and Multistage Sampling

In the theoretical formulation of simple random sampling and stratified random sampling, the statistical inquiry presumes the immediate availability of an exhaustive, pristine sampling frame listing every individual observation unit (such as a person, patient, or household). In large-scale national, regional, and industrial surveys, this ideal assumption fails completely. Constructing and continuously updating a comprehensive register of tens of millions of dispersed elements across an entire subcontinent is often financially impossible, logistically prohibitive, and outdated before field enumeration can even commence. Moreover, even if such an element-level frame were theoretically obtainable, a pure simple random sample of individual households scattered across a massive geographical territory would require interviewers to travel exorbitant distances to conduct single, isolated interviews, inflating per-unit data collection costs to unsustainable levels.

To overcome these structural limitations, survey sampling theory introduces cluster sampling. Instead of selecting individual observation units directly from an exhaustive master list, the population is partitioned into identifiable, mutually exclusive, and collectively exhaustive groups termed clusters or Primary Sampling Units (PSUs). A sampling frame of these clusters is constructed—an operational task that is vastly simpler and less expensive than listing all elements, because clusters correspond to existing administrative or physical entities such as census enumeration areas, municipalities, villages, school districts, or industrial plants. A probability sample of clusters is then selected.

In single-stage cluster sampling, every observation unit within each selected cluster is completely enumerated. In multistage sampling (or sub-sampling), the selection process proceeds hierarchically through successive tiers: first, a sample of primary sampling units is drawn; subsequently, within each selected primary unit, a probability sample of Secondary Sampling Units (SSUs) (such as dwellings or classrooms) is selected; in three-stage designs, ultimate Tertiary Sampling Units (TSUs) (such as individuals or students) are sampled within the secondary units.

Cluster sampling therefore resolves two paramount operational bottlenecks:

1. Frame Construction Efficiency: A complete list of all ultimate elements in the population is never required. The statistician requires only a listing of primary units at the initial stage. Detailed listings of secondary units are compiled exclusively within the small subset of primary units that are actually selected into the sample.

2. Fieldwork and Logistical Economy: By concentrating interviews within geographically compact clusters, travel time, fuel expenditures, lodging overhead, and field supervision logistics are reduced by orders of magnitude compared to unclustered, dispersed sampling designs.

However, this substantial operational efficiency comes at an inevitable statistical price. In human, ecological, and economic populations, elements residing within the same cluster tend to resemble one another more closely than elements chosen at random from different clusters. Households in the same village share the same agro-climatic conditions, local markets, public infrastructure, and socioeconomic status. Students in the same classroom share the same teacher, curriculum, and physical learning environment. Consequently, observing multiple elements from the same cluster provides partially redundant information. This internal homogeneity induces a positive correlation among units within clusters, leading to an inflation of sampling variance relative to simple random sampling of the same total number of elements. The central mathematical challenge of cluster sampling is to formalize, quantify, and minimize this variance inflation through rigorous design choices.

We begin by establishing the formal notation and mathematical framework for single-stage cluster sampling where all clusters contain an identical number of observation units. While equal cluster sizes represent an idealized scenario, this classical formulation provides the clearest mathematical foundation for understanding variance decomposition and intracluster correlation.

Population Structure and Notation

Consider a finite population $U$ partitioned into $N$ distinct clusters, with each cluster containing exactly $M$ observation units. The total number of elements in the entire finite population is therefore:

Let $y_{ij}$ denote the value of the survey characteristic for the $j$-th observation unit ($j = 1, 2, \dots, M$) within the $i$-th cluster ($i = 1, 2, \dots, N$).

We define the primary population parameters:

1. Cluster Total for Cluster $i$:

2. Cluster Mean for Cluster $i$:

3. Overall Population Total:

4. Overall Population Mean per Observation Unit:

Notice that because all clusters have identical size $M$, the overall population mean per element $\bar{\bar{Y}}$ is exactly equal to the simple unweighted average of the $N$ cluster means $\bar{Y}_i$.

Population Variance Decomposition

The total variability of the characteristic $y$ across all $NM$ elements in the finite population is quantified by the overall population variance:

In cluster sampling, this total variation is decomposed into two orthogonal components: variation between cluster means (between-cluster variance) and variation among elements within the same cluster (within-cluster variance).

To formalize this decomposition, we express the total sum of squares ($\text{SST}$) as:

Expanding the squared bracket:

The cross-product term evaluates to zero identically because $\sum_{j=1}^M (y_{ij} - \bar{Y}_i) = 0$ for every cluster $i$:

Therefore, we have the exact algebraic identity:

which is symbolically expressed as:

where $\text{SSB}$ is the between-cluster sum of squares and $\text{SSW}$ is the within-cluster sum of squares.

We define the corresponding mean squares:

1. Between-Cluster Mean Square (per element):

Note that the variance of the cluster totals $Y_i = M \bar{Y}_i$ is $S_t^2 = M^2 S_b^2$.

2. Within-Cluster Mean Square:

where $S_i^2 = \frac{1}{M - 1} \sum_{j=1}^M (y_{ij} - \bar{Y}_i)^2$ is the variance among elements within cluster $i$.

The relation between the overall variance $S^2$, between-cluster mean square $S_b^2$, and within-cluster mean square $S_w^2$ is obtained by dividing the sum of squares identity by $NM - 1$:

For large populations where $N$ and $NM$ are sufficiently large so that $N-1 \approx N$ and $NM-1 \approx NM$, dividing through by $NM$ yields the familiar asymptotic relationship:

The degree of internal similarity among units belonging to the same cluster is mathematically captured by the intracluster correlation coefficient (also known as the intraclass correlation coefficient), denoted by $\rho$. Understanding $\rho$ is essential to mastering cluster sampling, as it dictates the entire magnitude of design efficiency and variance inflation.

Formal Definition of $\rho$

In finite population sampling, the intracluster correlation coefficient $\rho$ is defined as the Pearson product-moment correlation coefficient between pairs of distinct elements belonging to the same cluster, evaluated across all possible distinct pairs in the entire population.

In cluster $i$, there are $M$ elements, yielding $M(M - 1)$ ordered pairs of distinct elements $(y_{ij}, y_{ik})$ with $j \ne k$. Summing across all $N$ clusters, there are exactly $N M (M - 1)$ such ordered pairs. The intracluster correlation coefficient is defined as:

To express $\rho$ in terms of the between-cluster and within-cluster mean squares, we manipulate the numerator. For any cluster $i$:

Recognizing that $\sum_{j=1}^M (y_{ij} - \bar{\bar{Y}}) = M (\bar{Y}_i - \bar{\bar{Y}})$, we have:

Summing both sides across all $N$ clusters:

Substituting our sum of squares definitions:

Thus, the numerator of $\rho$ is precisely:

Dividing this by the denominator $(M - 1)(NM - 1) S^2$, we obtain the exact algebraic formulation for $\rho$:

Rearranging this formula gives an exact expression for the between-cluster variance $S_b^2$ in terms of the total variance $S^2$ and the intracluster correlation $\rho$:

When $N$ is reasonably large such that $\frac{NM - 1}{M(N - 1)} \approx 1$, this simplifies to the fundamental relationship:

Mathematical Bounds of $\rho$

Unlike a standard bivariate correlation coefficient, which ranges over the symmetric interval $[-1, 1]$, the intracluster correlation coefficient $\rho$ is asymmetric:

To prove this fundamental result:

1. Upper Bound ($\rho \le 1$):

The between-cluster mean square $S_b^2$ is bounded above by the total variance scaled by $M$. If all units within every cluster are identical (complete internal homogeneity), then for every cluster $i$, $y_{ij} = \bar{Y}_i$ for all $j$. In this extreme scenario, $S_w^2 = 0$, so $\text{SSW} = 0$, which implies $\text{SST} = \text{SSB}$. Then:

Substituting this into the definition of $\rho$:

Hence, the maximum value of $\rho$ is $+1$, achieved when clusters are completely internally homogeneous.

2. Lower Bound ($\rho \ge -\frac{1}{M - 1}$):

The between-cluster mean square $S_b^2$ is a sum of squared deviations and must satisfy $S_b^2 \ge 0$. Setting $S_b^2 = 0$ corresponds to the case where all cluster means are identical: $\bar{Y}_1 = \bar{Y}_2 = \dots = \bar{Y}_N = \bar{\bar{Y}}$. Under this condition:

Therefore, the minimum possible value of $\rho$ is $-\frac{1}{M - 1}$.

Substantive Interpretation in Empirical Practice

In practical survey work, what values does $\rho$ assume?

• When $M$ is large (for example, census enumeration areas containing $M = 100$ or $200$ households), the lower bound $-\frac{1}{M - 1}$ is practically zero (e.g., $-0.01$ or $-0.005$). Thus, negative values of $\rho$ are virtually impossible in large clusters.
• Because human settlements, natural environments, and socioeconomic characteristics naturally cluster in space, elements within the same geographic cluster are systematically more similar to one another than to elements from distant clusters. Therefore, in empirical surveys, $\rho$ is almost universally positive.
• Even seemingly small values of $\rho$ (e.g., $\rho = 0.05$ or $0.10$) have devastating effects on precision when cluster size $M$ is substantial. For instance, if $M = 101$ and $\rho = 0.05$, the multiplier $1 + (M - 1)\rho = 1 + 100(0.05) = 6.0$. This means the sampling variance is inflated by a factor of 6 compared to simple random sampling.

Suppose a simple random sample without replacement (SRSWOR) of $n$ clusters is selected from the $N$ clusters in the population. Every unit within each selected cluster is enumerated, yielding a total observed sample of $n M$ elements.

The Sample Mean Estimator

Let $s \subset \{1, 2, \dots, N\}$ denote the set of $n$ sampled cluster indices. The sample mean per observation unit is:

where $\bar{y}_i = \frac{1}{M} \sum_{j=1}^M y_{ij} = \bar{Y}_i$ is the exact mean of the $i$-th cluster, since all units in the selected cluster are observed without error.

Because the clusters are selected via SRSWOR from $N$ clusters, $\bar{y}_{cl}$ is formally equivalent to the sample mean of $n$ observations drawn from a finite population of $N$ cluster means $\{\bar{Y}_1, \dots, \bar{Y}_N\}$.

Expectation and Design Unbiasedness

Taking the design expectation under the randomization distribution:

Thus, the cluster sample mean $\bar{y}_{cl}$ is an strictly design-unbiased estimator of the overall population mean per element $\bar{\bar{Y}}$.

Similarly, the estimator for the population total $Y$ is:

which is strictly design-unbiased: $E(\hat{Y}_{cl}) = Y$.

Sampling Variance Derivation

Applying the classical finite-population variance formula for SRSWOR of $n$ units from $N$ units:

where $f = \frac{n}{N}$ is the primary sampling fraction, and $S_b^2 = \frac{1}{N - 1} \sum_{i=1}^N (\bar{Y}_i - \bar{\bar{Y}})^2$ is the between-cluster mean square.

Substituting our relation between $S_b^2$, $S^2$, and $\rho$:

For large $N$, this expression simplifies to the fundamental equation:

Now compare this with the variance of a simple random sample without replacement of the same total number of elements, $m_{tot} = n M$:

where $f_{elem} = \frac{n M}{N M} = \frac{n}{N} = f$.

The ratio of the two variances is the Design Effect ($\text{Deff}$):

This celebrated formula, popularized by Leslie Kish (1965), proves that:

1. If $\rho = 0$ (elements within clusters are as heterogeneous as the general population), $\text{Deff} = 1$, and cluster sampling is just as efficient as simple random sampling.
2. If $\rho > 0$, $\text{Deff} > 1$, indicating variance inflation. The inflation grows linearly with the cluster size $M$.
3. If $\rho < 0$ (elements within clusters are more heterogeneous than the population at large—an artificial condition rarely found in geography), $\text{Deff} < 1$, and cluster sampling is actually more precise than simple random sampling.

Effective Sample Size

Because the design effect quantifies variance inflation, it reveals how many independent, uncorrelated observations the cluster sample is actually worth. The effective sample size ($n_{eff}$) is defined as:

For example, if a survey samples $n = 50$ clusters of size $M = 20$ (total elements $n M = 1,000$) and the intracluster correlation is $\rho = 0.10$:

Although the survey interviewers enumerate $1,000$ individuals, the statistical information gathered is equivalent to only $345$ individuals selected through pure simple random sampling. The remaining $655$ interviews represent informational redundancy caused by within-cluster correlation.

Unbiased Estimation of Variance from Sample Data

In practice, the true between-cluster mean square $S_b^2$ is unknown. We estimate it from the sample using the sample between-cluster mean square:

Since the $n$ cluster means are selected by SRSWOR from the $N$ cluster means, $s_b^2$ is an strictly unbiased estimator of $S_b^2$:

Consequently, the unbiased estimator of the variance of the cluster sample mean is:

Notice the remarkable simplicity of this estimator: to compute the standard error of a single-stage cluster sample, one does not need to compute within-cluster variances. One treats the observed cluster totals or cluster means as the fundamental data points and applies the standard SRS variance formula directly at the cluster level.

In real-world surveys, natural clusters (such as villages, city blocks, hospital wards, or schools) rarely contain equal numbers of observation units. Let $M_i$ denote the number of observation units in cluster $i$ ($i = 1, 2, \dots, N$). The total population size is:

and the average cluster size is $\bar{M} = \frac{M_0}{N}$.

When cluster sizes $M_i$ vary, several competing estimators of the population mean per element $\bar{\bar{Y}} = \frac{Y}{M_0}$ can be constructed, each with distinct statistical properties.

1. The Simple Mean of Cluster Means (Unweighted Estimator)

A naive estimator is the simple average of the sample cluster means:

where $\bar{y}_i = \frac{1}{M_i} \sum_{j=1}^{M_i} y_{ij} = \bar{Y}_i$.

The expectation of this estimator is:

Notice that $E(\bar{y}_u)$ is equal to $\bar{\bar{Y}} = \frac{\sum_{i=1}^N M_i \bar{Y}_i}{\sum_{i=1}^N M_i}$ if and only if the cluster size $M_i$ is completely uncorrelated with the cluster mean $\bar{Y}_i$. In economic and social applications, cluster size is frequently correlated with the attribute of interest (e.g., larger cities have higher average household incomes or higher crime rates). Therefore, $\bar{y}_u$ is generally biased for $\bar{\bar{Y}}$:

Unless the covariance between cluster size and cluster mean is known to be zero, this unweighted estimator should be avoided.

2. The Design-Unbiased Horvitz–Thompson Estimator

To obtain a strictly design-unbiased estimator when clusters are selected with equal probability (SRSWOR), we first estimate the population total $Y$ by scaling the sample cluster totals:

where $Y_i = \sum_{j=1}^{M_i} y_{ij}$ is the total for cluster $i$.

Dividing by the known total population size $M_0$:

Because $E\left[ \frac{N}{n} \sum_{i \in s} Y_i \right] = \sum_{i=1}^N Y_i = Y$, this estimator is strictly unbiased:

The exact variance of $\hat{\bar{Y}}_{unb}$ under SRSWOR of clusters is:

where $S_u^2 = \frac{1}{N - 1} \sum_{i=1}^N (Y_i - \bar{M} \bar{\bar{Y}})^2 = \frac{1}{N - 1} \sum_{i=1}^N (Y_i - \bar{Y}_t)^2$, and $\bar{Y}_t = \frac{Y}{N}$ is the average cluster total.

While strictly unbiased, $\hat{\bar{Y}}_{unb}$ suffers from a severe practical drawback: its variance depends on the variability of the cluster totals $Y_i$. Since $Y_i = M_i \bar{Y}_i$, wide variation in cluster sizes $M_i$ induces massive variation in $Y_i$, resulting in an unacceptably large sampling variance. Furthermore, calculating this estimator requires knowing the exact population total size $M_0$, which is frequently unknown when conducting surveys in developing or dynamic settings.

3. The Ratio-to-Size Estimator (The Combined Cluster Mean)

The standard, most widely adopted estimator in survey practice is the ratio estimator:

This estimator divides the total sample aggregate of the characteristic $y$ by the total number of elements actually observed in the sample clusters. It is a ratio of two random variables: the estimated total of $y$ and the estimated total number of elements.

Key Advantages of $\hat{\bar{Y}}_R$:

1. Self-Weighting and Operational Feasibility: It does not require knowledge of the total population size $M_0$.

2. Substantial Variance Reduction: The cluster total $Y_i$ is typically strongly correlated with cluster size $M_i$. By dividing by $\sum_{i \in s} M_i$, the ratio estimator absorbs the fluctuations in cluster sizes, dramatically dampening sampling variability.

Statistical Properties:

Like all ratio estimators, $\hat{\bar{Y}}_R$ is technically biased in finite samples, but its bias is of order $O(1/n)$ and becomes negligible when the number of sampled clusters $n$ exceeds 30.

Using a first-order Taylor series approximation, the mean squared error (approximate variance) is:

An approximately unbiased sample estimator of this variance is:

where $\bar{m} = \frac{1}{n} \sum_{i \in s} M_i$ is the average cluster size observed in the sample.

When cluster sizes $M_i$ vary drastically, simple random sampling of clusters followed by ratio estimation is still sub-optimal. The modern standard methodology is to select clusters with Probability Proportional to Size (PPS), which we investigate thoroughly in Unit 5.

When primary sampling units are large (such as entire districts, cities, or large health facility catchment areas), surveying every single element within selected PSUs is neither affordable nor necessary. Instead, the statistician implements two-stage cluster sampling:

• Stage 1: Select a probability sample of $n$ primary sampling units (PSUs) out of $N$.
• Stage 2: Within each selected PSU $i \in s$, select a probability sample of $m_i$ secondary sampling units (SSUs) out of the $M_i$ available elements.

Conditional Expectation and the Law of Total Variance

The mathematical analysis of multistage sampling relies on the fundamental laws of iterated conditional expectation and conditional variance. Let $\hat{\theta}$ be an estimator of a population parameter $\theta$, and let $s_1$ denote the specific sample of PSUs chosen at the first stage.

1. Law of Iterated Expectation:

where $E_2(\cdot \mid s_1)$ denotes expectation over all possible second-stage samples conditional on the fixed first-stage sample $s_1$, and $E_1(\cdot)$ denotes expectation over all possible first-stage samples.

2. Law of Total Variance (Variance Decomposition Theorem):

This decomposition partitions the total variance into two distinct structural mechanisms:

• First-Stage Variance Component ($V_1$): The variability arising because different sets of PSUs could have been chosen at the first stage.
• Second-Stage Variance Component ($E_1 V_2$): The average variability introduced by sampling a subset of SSUs within the selected PSUs, rather than enumerating them completely.

Two-Stage Sampling with Equal Cluster and Sub-Sample Sizes

Consider a population of $N$ clusters, each containing $M$ elements. At the first stage, $n$ clusters are selected using SRSWOR ($f_1 = n/N$). At the second stage, $m$ elements are selected from each chosen cluster using SRSWOR ($f_2 = m/M$). The overall sampling fraction is $f = f_1 f_2 = \frac{n m}{N M}$.

The overall sample mean is:

where $\bar{y}_i = \frac{1}{m} \sum_{j=1}^m y_{ij}$ is the sample mean of the $m$ elements drawn from cluster $i$.

Unbiasedness:

First, condition on the selected primary units $s_1$:

Taking the first-stage expectation over all possible selections of $s_1$:

Thus, $\bar{\bar{y}}$ is strictly design-unbiased.

Derivation of the Exact Sampling Variance:

Using the Law of Total Variance:

1. First-Stage Component:

Since $E_2(\bar{\bar{y}} \mid s_1) = \frac{1}{n} \sum_{i=1}^n \bar{Y}_i$, this is the sample mean of $n$ units drawn by SRSWOR from the population of $N$ cluster means. Therefore:

where $S_b^2 = \frac{1}{N - 1} \sum_{i=1}^N (\bar{Y}_i - \bar{\bar{Y}})^2$.

2. Second-Stage Component:

Conditional on $s_1$, the sub-samples drawn within different PSUs are completely independent. Therefore, the variance of the sum is the sum of the variances:

Within cluster $i$, an SRSWOR of $m$ elements is drawn from $M$ elements. Hence:

where $S_i^2 = \frac{1}{M - 1} \sum_{j=1}^M (y_{ij} - \bar{Y}_i)^2$. Substituting this into the sum:

Now, take the expectation $E_1$ over all first-stage samples. Since the $n$ clusters are chosen by SRSWOR from $N$:

Therefore:

Combining both components yields the celebrated Two-Stage Variance Formula:

Structural Analysis of the Variance Formula:

Notice how the formula behaves under limiting design choices:

• If $m = M$, then $f_2 = 1$, and the second term vanishes completely:

This matches our earlier single-stage cluster sampling variance exactly.

• If $n = N$, then $f_1 = 1$, and the first term vanishes:

This matches stratified random sampling with proportional allocation, where every cluster serves as a stratum.

• In two-stage cluster sampling, both terms are present. The first term represents the penalty for not sampling all clusters; the second term represents the penalty for not enumerating all elements within the sampled clusters.

Sample Variance Estimation in Two-Stage Sampling

To estimate $V(\bar{\bar{y}})$ from empirical sample data, we define two sample mean squares:

1. Sample Between-Cluster Mean Square:

2. Sample Within-Cluster Mean Square:

where $s_i^2 = \frac{1}{m - 1} \sum_{j=1}^m (y_{ij} - \bar{y}_i)^2$.

It is straightforward to show that $s_2^2$ is an unbiased estimator of the population within-cluster variance:

However, $s_1^2$ is not an unbiased estimator of $S_b^2$. Because $\bar{y}_i$ contains second-stage sampling error as an estimate of $\bar{Y}_i$, $s_1^2$ is inflated:

Solving for an unbiased estimator of $V(\bar{\bar{y}})$:

The Ultimate Cluster Principle

Notice what happens when the first-stage sampling fraction $f_1 = \frac{n}{N}$ is small (for example, $f_1 < 0.05$, which is typical in large-scale national surveys where 300 enumeration areas are selected out of 50,000 national enumeration areas): The term $\frac{f_1(1 - f_2)}{nm} s_2^2$ becomes negligibly small, and $1 - f_1 \approx 1$. The variance estimator simplifies to:

This foundational result is known as the Ultimate Cluster Principle (Hansen, Hurwitz, and Madow, 1953). When the first-stage sampling fraction is small, the complex internal multistage sampling structure within PSUs can be completely ignored when estimating total variance! The variation among the PSU totals or PSU means automatically captures both the first-stage and second-stage variance components.

A central operational question in designing a two-stage survey is: Given a fixed total budget, what combination of primary sampling units $n$ and secondary units per PSU $m$ will minimize the sampling variance? Or conversely, to achieve a pre-specified target variance, what allocation minimizes total financial cost?

Linear Cost Model

In field operations, survey costs are accurately characterized by a linear cost function:

where:

• $c_0$ is the fixed overhead cost (planning, frame development, software, administrative management).
• $c_1$ is the cost per Primary Sampling Unit (traveling to the cluster, locating boundaries, listing households, setting up field accommodation).
• $c_2$ is the marginal cost per Secondary Sampling Unit (conducting the interview with a household, administering questionnaires, biometric measurement).
• $C - c_0 = C_{field} = n(c_1 + m c_2)$ is the net operational field budget.

Optimal Allocation Formula

Ignoring finite population corrections ($f_1 \approx 0, f_2 \approx 0$), the variance of the two-stage mean is:

and the net field cost is:

We seek to minimize $V$ subject to fixed $C_{field}$, or minimize $C_{field}$ subject to fixed $V$. By the Cauchy–Schwarz inequality, the product of variance and cost is minimized when:

Differentiating with respect to $m$ and setting to zero:

Taking the positive square root yields the Optimal Sub-Sample Size:

Using our relationship with the intracluster correlation $\rho$, where $\frac{S_w^2}{S_b^2} \approx \frac{1 - \rho}{\rho}$:

Once $m_{opt}$ is determined, the optimal number of clusters $n_{opt}$ is determined by the budget constraint:

Practical Interpretation of $m_{opt}$

The mathematical structure of $m_{opt}$ reveals two fundamental design rules:

1. Cost Ratio Effect ($\sqrt{c_1 / c_2}$): If visiting a cluster is very expensive relative to conducting an individual interview ($c_1 \gg c_2$), $m_{opt}$ increases. It pays to take more interviews per cluster to amortize the steep initial travel cost.

2. Correlation Effect ($\sqrt{(1 - \rho) / \rho}$): If the intracluster correlation $\rho$ is high, $m_{opt}$ decreases. Because units within the cluster are highly redundant, sampling many units within the same cluster wastes resources; the budget is better spent visiting more clusters ($n \uparrow$) and taking few interviews per cluster ($m \downarrow$).

Systematic sampling is one of the most widely applied sampling techniques in survey practice. In systematic sampling, elements are selected at equal numeric intervals from an ordered list.

Selection Protocol and Notation

Suppose a population of $N$ units is indexed from $1$ to $N$. We wish to draw a sample of $n$ units. Let $k$ be an integer such that:

The integer $k$ is called the sampling interval.

The selection protocol proceeds as follows:

1. Draw a single random integer $r$ uniformly from the set $\{1, 2, \dots, k\}$. The value $r$ is the random start.
2. The sample consists of the units with indices:

Equivalence to Single-Stage Cluster Sampling

Statistically, systematic sampling is mathematically identical to selecting exactly one cluster out of $k$ available clusters with equal probability $1/k$.

To see this correspondence, consider arranging the $N$ units of the population in a rectangular matrix of $k$ rows and $n$ columns:

Each row in this matrix represents a potential systematic sample. There are exactly $k$ possible samples, each containing $n$ elements. Selecting a random start $r \in \{1, \dots, k\}$ is exactly equivalent to selecting one row (one cluster) at random from a population of $k$ clusters, where each cluster contains $n$ elements!

Variance of the Systematic Sample Mean

Let $\bar{y}_{sys} = \bar{y}_{r\cdot} = \frac{1}{n} \sum_{j=0}^{n-1} y_{r + jk}$ denote the sample mean of the systematic sample.

Because each of the $k$ rows is selected with equal probability $1/k$:

Thus, the systematic sample mean is strictly design-unbiased.

The variance of $\bar{y}_{sys}$ over all $k$ possible samples is:

Using our earlier intracluster correlation identity for clusters of size $n$, we can express the variance of $\bar{y}_{sys}$ in terms of the population variance $S^2$ and the intracluster correlation among elements within the same systematic sample, denoted by $\rho_{sys}$:

where $\rho_{sys}$ is the correlation between pairs of elements in the same systematic sample:

Comparison of Systematic Sampling with SRS and Stratified Sampling

The performance of systematic sampling depends critically on how the population list is ordered:

1. Randomly Ordered Population: If the frame elements are ordered randomly (e.g., alphabetically by an unrelated name), $\rho_{sys} \approx -\frac{1}{N - 1} \approx 0$. Then $V(\bar{y}_{sys}) \approx V_{srs}(\bar{y})$. Systematic sampling performs identically to SRSWOR, but is operationally simpler to execute.

2. Monotonic Trend (Linear Trend): If the list is ordered monotonically with respect to the survey variable (e.g., sorting farms by acreage or employees by salary), each systematic sample contains elements from the bottom, middle, and top of the distribution. In this case, within-sample heterogeneity is maximized, forcing $\rho_{sys} < 0$. Consequently:

Systematic sampling behaves like an implicitly stratified sample with one unit per stratum of size $k$, achieving outstanding precision.

3. Periodic Variation (The Fatal Trap): If the population exhibits cyclical or periodic patterns with a period $P$ that coincides with the sampling interval $k$ (or a multiple thereof), every selected element lands on the exact same phase of the cycle (e.g., always sampling retail sales on Mondays, or selecting corner lots in a grid city). In this catastrophic scenario, elements within the systematic sample are nearly identical, driving $\rho_{sys} \to 1$. The variance explodes, and $\bar{y}_{sys}$ becomes severely inefficient.

In Section 4.4, we derived the two-stage cluster sampling variance under the assumption that all $N$ primary clusters contained an identical number of observation units $M$, and that an equal number of elements $m$ was sampled from each selected cluster. In real-world demographic, agricultural, and socioeconomic surveys, clusters are natural administrative or geographical entities—villages, urban wards, schools, census enumeration areas—whose population sizes $M_i$ vary by orders of magnitude.

When cluster sizes $M_i$ vary, the statistical design must formalize both the first-stage selection of primary sampling units (PSUs) and the second-stage selection of secondary sampling units (SSUs) to maintain strict design consistency and avoid uncontrolled variance inflation.

The General Two-Stage Horvitz–Thompson Estimator

Consider a finite population $U$ partitioned into $N$ primary sampling units $U_1, U_2, \dots, U_N$. PSU $i$ contains $M_i$ elements. The total number of elements in the entire finite population is:

Let $y_{ij}$ denote the survey measurement for element $j \in U_i$ ($j = 1, 2, \dots, M_i$). The total for PSU $i$ is $Y_i = \sum_{j=1}^{M_i} y_{ij}$, and the overall population total is $Y = \sum_{i=1}^N Y_i$.

The two-stage sampling protocol executes as follows:

• Stage 1: Select a sample $s_1$ of $n$ PSUs from $\{1, 2, \dots, N\}$ according to a probability design without replacement characterized by first-order inclusion probabilities $\pi_i$ and second-order joint inclusion probabilities $\pi_{ik}$ ($i, k = 1, \dots, N$).
• Stage 2: Within each selected PSU $i \in s_1$, select an SRSWOR sample $s_{2i}$ of $m_i$ elements from the $M_i$ available units, with conditional inclusion probabilities $\pi_{j|i} = m_i / M_i$.

For each selected PSU $i \in s_1$, the design-unbiased estimator of the PSU total $Y_i$ is:

where $\bar{y}_i = \frac{1}{m_i} \sum_{j \in s_{2i}} y_{ij}$ is the sample mean of the elements observed in PSU $i$.

The Two-Stage Horvitz–Thompson Estimator of the overall population total $Y$ is:

To estimate the overall population mean per element $\bar{\bar{Y}} = Y / M_0$:

Proof of Strict Design Unbiasedness

Using the Law of Iterated Expectation:

Conditioning on the selected primary units $s_1$:

Taking the first-stage expectation over all possible selections of $s_1$:

This proves that the two-stage Horvitz–Thompson estimator is strictly design-unbiased for any cluster sizes $M_i$ and any sub-sample sizes $m_i$.

Exact Variance Decomposition for Unequal Cluster Sizes

Applying the Law of Total Variance:

1. First-Stage Variance Component ($V_1$):

Since $E_2(\hat{Y}_{2HT} \mid s_1) = \sum_{i \in s_1} \frac{Y_i}{\pi_i}$, this is simply the Horvitz–Thompson estimator of the total of the population of cluster totals $\{Y_1, \dots, Y_N\}$. Under the Sen–Yates–Grundy formulation (for fixed $n$):

2. Second-Stage Variance Component ($E_1 V_2$):

Conditional on $s_1$, sub-sampling across different PSUs is completely independent:

Within PSU $i$, an SRSWOR of $m_i$ elements is drawn from $M_i$ elements. The variance of the estimated total $\hat{Y}_i = M_i \bar{y}_i$ is:

where $f_{2i} = m_i / M_i$ and $S_i^2 = \frac{1}{M_i - 1} \sum_{j=1}^{M_i} (y_{ij} - \bar{Y}_i)^2$. Substituting this into the sum:

Taking the first-stage expectation $E_1$:

Combining both components yields the general Two-Stage Variance Formula for Unequal Clusters:

Unbiased Sample Variance Estimation

To estimate $V(\hat{Y}_{2HT})$ from empirical sample data, we apply the Sen–Yates–Grundy estimator to the first stage, substituting the estimated cluster totals $\hat{Y}_i$ for the true totals $Y_i$, and subtracting the second-stage noise that inflates the between-cluster mean square.

The strictly unbiased sample variance estimator is:

where $s_i^2 = \frac{1}{m_i - 1} \sum_{j \in s_{2i}} (y_{ij} - \bar{y}_i)^2$ is the sample variance of the elements observed in PSU $i$.

When the first-stage sampling fraction is small ($\pi_i \ll 1$), the first term dominates, and the Ultimate Cluster estimator provides an extraordinarily accurate, simple approximation.

In our initial investigation of optimal sample allocation (Section 4.5), survey operational costs were characterized by a linear model: $C = c_0 + n c_1 + n m c_2$. Under this simple model, the marginal cost of visiting a PSU ($c_1$) was assumed to be independent of the number of PSUs selected ($n$).

In actual field operations across large geographic expanses (such as provinces, agricultural river basins, or forestry tracts), the linear travel cost assumption fails. When field interview teams travel between sampled PSUs, the distance traveled between clusters depends directly on how densely or sparsely the $n$ clusters are scattered across the geographic territory.

The Beardwood–Halton–Hammersley Spatial Travel Law

In mathematical spatial optimization, the Traveling Salesman Problem (TSP) addresses the shortest closed tour connecting $n$ points distributed across a two-dimensional geographic region of surface area $A$.

A landmark theorem in stochastic geometry proved by J. Beardwood, J. H. Halton, and J. M. Hammersley (1959) establishes that if $n$ points are distributed uniformly at random across a bounded region of area $A$, the minimum travel distance $D_n$ required to visit all $n$ locations satisfies:

where $\beta_{TSP}$ is the universal Traveling Salesman constant for Euclidean two-dimensional space.

Consequently, the total mileage or travel distance between clusters does not grow linearly with $n$, but scales proportionally to the square root of the number of clusters:

The Non-Linear Spatial Survey Cost Model

Incorporating the Beardwood–Halton–Hammersley spatial law into survey operations yields the Non-Linear Spatial Cost Function (Hansen, Hurwitz, and Madow, 1953; Groves, 1989):

where:

• $c_0$ is fixed administrative overhead.
• $c_t$ is the travel cost per unit distance (fuel expenditures, vehicle maintenance, travel time wages).
• $A$ is the total geographical area of the survey universe (in square kilometers).
• $c_1$ is the local fixed cost per PSU (locating boundary monuments, meeting local community leaders, establishing field headquarters).
• $c_2$ is the marginal interview cost per secondary unit (household interview time, respondent incentives).

The term $c_t \sqrt{n A}$ represents the inter-cluster travel cost, while $n c_1 + n m c_2$ represents the intra-cluster operational cost.

Optimal Allocation under Spatial Travel Constraints

To minimize the two-stage sampling variance $V \approx \frac{S_b^2}{n} + \frac{S_w^2}{n m}$ subject to the non-linear spatial budget constraint:

We set up the Lagrangian:

Differentiating with respect to $m$ and $n$:

1. With respect to $m$:

1. With respect to $n$:

Equating the expressions for $\lambda$ and solving for $m$:

Taking the positive square root yields the Spatially Adjusted Optimal Sub-Sample Size:

Methodological Insights from the Spatial Formulation:

1. Density Dampening: The effective cost per PSU is inflated by the spatial travel term $\frac{c_t \sqrt{A}}{2 \sqrt{n}}$. In vast, sparsely populated geographic regions (large area $A$), traveling between clusters is the dominant operational hurdle. This inflates the effective cost per PSU, forcing $m_{opt}$ upward: interview teams must conduct more interviews per cluster ($m \uparrow$) to amortize the grueling travel required to reach each isolated site.

2. Urban Compactness: In compact urban environments (small area $A$), inter-cluster travel is rapid and inexpensive ($c_t \sqrt{A} \to 0$). The optimal design shifts toward sampling many compact clusters with very few interviews per cluster ($n \uparrow, m \downarrow$), which minimizes variance inflation by suppressing the intracluster correlation penalty $1 + (m - 1)\rho$.

In Section 4.6, we proved that systematic sampling is formally identical to selecting a single cluster of size $n$ from $k$ systematic clusters. Because only one systematic cluster is observed, the classical textbook states that an unbiased estimator of variance does not exist without making strong modeling assumptions.

In modern survey practice, however, national statistical agencies frequently deploy Systematic Selection of Clusters (e.g., selecting census enumeration areas systematically along a geographical serpentining path). To provide rigorous, robust standard errors for systematic designs, Kirk Wolter (2007) formalized two primary replication frameworks:

1. The Successive Difference Replication Estimator

When clusters are ordered geographically or monotonically prior to systematic selection, adjacent clusters in the sample list are physically proximate and share similar environmental and demographic characteristics.

Let $\bar{y}_1, \bar{y}_2, \dots, \bar{y}_n$ denote the sequence of sample cluster means arranged in the exact order of their systematic selection. The Successive Difference Variance Estimator is defined as:

Proof of Trend Elimination:

Suppose the population values contain a deterministic linear trend across the ordered frame: $\bar{Y}_k = \alpha + \beta k$. The standard sample variance $s_b^2 = \frac{1}{n - 1} \sum (\bar{y}_k - \bar{y})^2$ is heavily inflated by the trend, since $\sum (\bar{Y}_k - \bar{Y})^2 \approx \frac{\beta^2 n^3}{12}$. Now evaluate the first differences:

The linear trend is completely transformed into a constant intercept $\beta$, which is subtracted out in higher-order differences! The Second-Order Successive Difference Estimator removes quadratic curvature:

Successive difference estimators protect survey standard errors against false inflation caused by spatial trends and geographic sorting.

2. Random Group Systematic Sampling

To achieve strictly unbiased variance estimation without sacrificing the operational simplicity of systematic sampling, Wolter recommended Random Group Systematic Selection:

1. Instead of selecting a single systematic sample with interval $k = N / n$, the desired total sample size $n$ is divided into $R$ independent sub-samples ($R \ge 5$ or $10$), each of size $m = n / R$.
2. The sampling interval is expanded to:

1. The statistician draws $R$ independent random starts $r_1, r_2, \dots, r_R$ from $\{1, 2, \dots, k'\}$ without replacement.
2. Each random start generates an independent systematic sample $s_r$ of size $m$, yielding $R$ independent systematic sample means $\bar{y}_1^*, \bar{y}_2^*, \dots, \bar{y}_R^*$.

The overall sample mean is the simple average:

The Random Group Variance Estimator is:

Because the $R$ random starts are drawn independently, $\hat{V}_{RG}$ is a strictly design-unbiased estimator of the true sampling variance! Random group systematic sampling completely resolves the non-estimability of systematic variance while preserving all the fieldwork efficiencies of systematic selection.

In designing multistage cluster surveys, one of the most critical decisions confronting the survey methodologist is determining the optimal cluster take (the number of households or individuals $m$ to interview within each sampled primary sampling unit). If $m$ is chosen too small (e.g. $m = 2$ or $3$), field travel costs explode because enumerator teams must journey across hundreds of geographically separated villages. Conversely, if $m$ is chosen too large (e.g. $m = 80$ or $100$), the survey suffers from severe variance inflation due to intracluster correlation, causing the design effect $\text{Deff} = 1 + (m - 1)\rho$ to skyrocket and wasting survey budget on redundant interviews.

To formalize this fundamental trade-off, Seymour Sudman (1976) and Hansen, Hurwitz, and Madow (1953) formulated mathematical models for the optimal cluster take $m^*$.

Mathematical Derivation of Sudman's Formula

Consider a two-stage cluster survey with equal cluster sizes $M$. The statistician wishes to estimate a population mean $\bar{Y}$. The total sampling variance under two-stage sampling is:

Now consider the classical field budget cost model:

where:

• $c_0$ is overhead administrative cost.
• $c_1$ is the cost of adding an additional primary sampling unit (cluster listing, community entry, supervisor travel, and lodging).
• $c_2$ is the marginal cost of interviewing an additional household within an already visited cluster.
• The total variable budget available for field operations is $C_{var} = C - c_0 = n (c_1 + c_2 m)$.

Solving the budget constraint for the number of clusters $n$:

Substitute this expression for $n$ into the variance formula:

Notice that $c_1 \rho$ and $c_2(1 - \rho)$ are constants independent of $m$. To find the value of $m$ that minimizes the variance $V(\bar{y})$ for a fixed budget $C_{var}$, differentiate with respect to $m$ and set to zero:

Equating the two terms:

Taking the positive square root yields Sudman's Optimal Cluster Take Formula:

Methodological Interpretation of Sudman's Rule

Sudman's formula reveals the exact mathematical interplay between cost structure and demographic clustering:

1. Ratio of Unit Costs ($c_1 / c_2$): The term $c_1 / c_2$ represents the relative cost of accessing a cluster compared to interviewing an individual within that cluster.

• If travel is cheap ($c_1 \approx c_2$), $m^*$ is small: the survey should visit many clusters and conduct very few interviews per cluster.
• If travel is extremely expensive (e.g., remote rural provinces where reaching a cluster requires four-wheel drive vehicles or boats, such that $c_1 / c_2 = 50$ or $100$), $m^*$ expands substantially, justifying larger cluster takes to amortize heavy travel costs.

2. Intraclass Correlation ($\rho$): The term $(1 - \rho)/\rho$ represents the demographic diversity of the cluster.

• When $\rho$ is small ($\rho \approx 0.01$), units within a cluster are nearly independent: $(1 - \rho)/\rho \approx 99$, driving $m^*$ upward.
• When $\rho$ is large ($\rho \approx 0.20$), units within a cluster are nearly identical copies of one another: $(1 - \rho)/\rho \approx 4$, forcing $m^*$ down to avoid collecting redundant data.

Numerical Walkthrough: Sudman's Rule in Field Practice

Suppose a national statistical agency is designing a household expenditure survey across rural districts:

• Cost to establish a cluster (cartography, listing, community leader consent, travel): $c_1 = 450 \text{ dollars}$.
• Cost per household interview: $c_2 = 18 \text{ dollars}$.
• Cost ratio: $c_1 / c_2 = 450 / 18 = 25.0$.
• Pilot survey estimate of intracluster correlation for household consumption expenditure: $\rho = 0.06$.

Evaluate Sudman's formula:

The optimal cluster take is exactly $20$ households per primary sampling unit. If the total available variable field budget is $C_{var} = 72,000 \text{ dollars}$, the optimal number of clusters to sample is:

The survey will visit $n = 89$ clusters, sampling $m = 20$ households per cluster, yielding a total sample size of $89 \times 20 = 1,780$ households. This allocation mathematically guarantees the lowest possible sampling variance for the fixed 72,000-dollar budget.

The practical execution of cluster and multistage surveys in developing and transitional statistical systems requires rigorous field management protocols to prevent non-sampling errors from overwhelming the theoretical gains of probability sampling.

Household Listing Protocols and Segmentation

When Primary Sampling Units are selected, field teams must enumerate the secondary sampling units. In rapidly growing urban centers or informal settlements, census cartography rapidly becomes obsolete. Survey agencies implement standardized listing protocols:

1. Canvassing and Boundary Verification: Field supervisors locate the physical boundaries of the selected PSU using GPS coordinates and satellite basemaps.

2. Segmentation: If the PSU contains more than 300 households, it is partitioned into 2 to 4 compact segments with well-defined physical boundaries, and one segment is selected with probability proportional to estimated size.

3. Numbering and Listing: Enumerators walk sequentially through the segment, assigning a unique identification number to every structure and dwelling unit, recording the name of the household head and physical address.

4. Equal Probability Systematic Selection: Once the complete listing of $M_i$ households is compiled, field teams apply a systematic sampling interval $k_i = M_i / m_i$ with a random start provided by the central office. Interviewers are strictly prohibited from selecting households on the spot, as informal substitution introduces severe selection bias.

Design Effect Monitoring and Weight Trimming

During survey analysis, survey statisticians calculate the empirical design effect for all primary survey indicators:

Key empirical guidelines include:

• Demographic indicators (e.g., age, sex ratio) typically exhibit low clustering: $\text{Deff} \approx 1.2 - 1.8$ ($\rho \approx 0.01 - 0.03$).
• Health and immunization indicators (e.g., clinic access, measles vaccination) exhibit moderate clustering: $\text{Deff} \approx 2.0 - 3.5$ ($\rho \approx 0.05 - 0.10$).
• Poverty, agricultural practices, and environmental conditions exhibit severe clustering: $\text{Deff} \approx 4.0 - 8.0+$ ($\rho \approx 0.15 - 0.30$).

When cluster sizes $M_i$ vary widely, sampling weights $w_{ij} = \frac{1}{\pi_{ij}}$ fluctuate substantially. Extreme weights inflate variance according to Kish's weight design effect formula:

Total design effect is approximately the product of clustering and weighting effects:

To protect survey precision against explosive variances caused by outlier weights, statistical agencies apply weight trimming, capping weights at the 95th or 99th percentile and recalibrating totals to preserve design consistency. In the next unit, we explore how unequal probability selection mechanisms—most notably Probability Proportional to Size (PPS)—eliminate this weight variance altogether.