
==== Front
MethodsX
MethodsX
MethodsX
2215-0161
Elsevier

S2215-0161(24)00379-0
10.1016/j.mex.2024.102928
102928
Statistic
On the implementation of stratified two-stage simple random sampling without replacement, with possible collapsed strata
Aubry Philippe philippe.aubry@ofb.gouv.fr

OFB - Office français de la biodiversité - Direction surveillance, évaluation, données - Unité données et appui méthodologique, Saint Benoist, BP 20, F-78612 Le Perray-en-Yvelines, France
03 9 2024
12 2024
03 9 2024
13 1029281 7 2024
10 8 2024
20 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Graphical abstract

Two-stage stratified sampling is a complex design that involves nested sampling units and stratification. This complexity increases when the strata have too few sampled units for variance estimation, necessitating the use of collapsed strata, where multiple strata are combined to ensure an adequate sample size. When collapsing strata, two cases can be distinguished depending on whether a size variable associated with the variable of interest is available at the stratum level.• We present computer-implementable formulas for total, mean, and ratio estimators, along with their corresponding sampling variance estimators, for stratified two-stage simple random sampling without replacement, and we provide ready-to-use algorithms.

• We introduce two methods for grouping strata: (1) a deterministic approach that uses stratum codes to define an ordinal variable, which orders the strata, and (2) a stochastic method that aims to minimize within-group inertia, which measures the heterogeneity within the newly formed groups of strata.

• We emphasize that, unlike the correlation between a size variable and the variable of interest at the stratum level, the bias of the sampling variance estimator for the collapsed strata technique is not invariant to linear transformations. It follows that a high correlation does not ensure a low-bias estimator of the sampling variance.

Keywords

Stratified two-stage sampling
Expansion estimator
Woodruff’s method
Collapsed strata
Size variable
Grouping strata
Combinatorial optimization
Method name

Computation in stratified two-stage simple random sampling without replacement, with possible collapsed strata
==== Body
pmcSpecifications tableSubject Area:	Computational statistics	
More specific subject area:	Complex sampling designs; Variance estimation techniques	
Method name:	Computation in stratified two-stage simple random sampling without replacement, with possible collapsed strata	
Name and reference of original method:	• M.H. Hansen, W.N. Hurwitz, W.G. Madow, Sample survey method and theory. Volume II. Theory. John Wiley & Sons, New York, USA, 1953.	
	• W.G. Cochran, Sampling techniques. Third edition. John Wiley & Sons, New York, USA, 1977.	
	• K. M. Wolter, Introduction to variance estimation. Second edition. Springer, New York, USA, 2007.	
Resource availability:	The example dataset in supplementary material can be used to reproduce the results presented in Section “Method validation”.	

Background

Stratified multistage sampling is widely recognized as the most common sampling design for large-scale surveys in a variety of fields, including socioeconomics, agriculture, natural resource assessment, and environmental monitoring.

Among stratified multistage sampling designs without replacement, the simplest is stratified two-stage simple random sampling (see Fig. 1). Although the theory behind this design is well established, it is often not presented in a form that is readily applicable in computer implementation. In addition, it is sometimes necessary to use the collapsed strata technique when the sample sizes within strata are too small to apply standard variance estimators or to obtain sufficiently precise estimates of within-strata variances. Similarly, although this technique is not new, it is often not presented in a way that is conducive to direct computer implementation. In particular, methods for grouping strata are rarely described or discussed in sufficient detail.Fig. 1 Schematic representation of stratified two-stage simple random sampling. (a) Eight primary sampling units (rectangles) are stratified into two strata (bright and light green backgrounds). Each primary sampling unit consists of a variable number of secondary sampling units (small circles). (b) In the first stage, primary sampling units are selected in each stratum by simple random sampling without replacement. (c) In the second stage, three secondary sampling units are selected in each primary sampling unit of the first-stage sample, again by simple random sampling without replacement.

Fig. 1

The first objective of this article is to transfer elements for implementing stratified two-stage sampling with simple random sampling and without replacement at both stages in terms of (i) the sampling algorithm; (ii) estimators for a total, a mean or a ratio; and (iii) (nonsimplified) estimators of their respective sampling variances. To estimate the sampling variance for a ratio, we rely on the estimators provided for the total, using the method of Woodruff [1]. Our perspective is operational rather than theoretical, and we provide the estimators to the reader without proof. Their algorithmic counterparts are given for ease of computer implementation in any procedural programming language. In addition, the reader is cautioned against possible misinterpretations about the nature of the variance estimators. For theoretical considerations, we refer the reader to [2, Ch. 6] or [3, Ch. 11], and for a general and modern presentation, we refer the reader to [4, Ch. 4] or [5, Ch. 7].

Similarly, for the collapsed strata technique, the goal of this article is to provide formulas with notations tailored for computer implementation. A significant part of the article is devoted to the introduction of two methods for grouping strata on the basis of available information. One method is deterministic, whereas the other takes a stochastic approach for the combinatorial optimization of a given objective function. We provide ready-to-use algorithms so that their computer implementation can be easily reproduced by readers in the programming language of their choice. For the theory of collapsed strata estimation, we refer the reader to [2,Ch. 9, Sec. 5], [3, Sec. 5A.12] and [6, Sec. 2.5].

This technical article is intended to provide a compendium of formal expressions and algorithms that may be useful to researchers or engineers who face similar needs and aims to help them save time in their own work. As a complement to Aubry [7], this article also has potential pedagogical value for graduate courses in computational statistics, contributing to in-depth insights into multistage sampling, given the ubiquity of stratified two-stage sampling (with possible collapsed strata) in large-scale surveys (see Aubry et al. [8] for a recent application). While our article does not aim to present new advances in survey sampling theory, its purpose is also to shed light on several overlooked aspects that we believe are important.

Method details

Stratified two-stage simple random sampling

The stratified two-stage sampling design combines a two-stage sampling design with a stratified sampling design in the first stage. This design implements a three-level hierarchy from the finite universe partitioned by the strata to the elementary or ultimate sampling units, which are sampling units that are not further disaggregated, in terms of which measurements or observations are made. These three levels are (i) a first level composed of strata, i.e., groups of homogeneous units with respect to a criterion to be defined (the strata are not sampled; they are all represented in the final sample selected); (ii) a second level composed of primary sampling units (PSUs), which group together units of the lower level; and (iii) a third level composed of secondary sampling units (SSUs), which are the elementary sampling units included in the PSUs (Fig. 1).

In this article, we use simple random sampling without replacement (SRSWOR) to sample the PSUs in the strata and the SSUs in the PSUs. Thus, in each stratum, a sample of PSUs is selected via SRSWOR; this is the first stage of sampling. In each PSU selected in the first stage, a sample of SSUs is selected via SRSWOR; this is the second stage of sampling (Fig. 1). Sampling in one stratum is independent of sampling in another stratum, and the same is true for sampling in the PSUs (independence; see, for example, [4, pp. 101-102, 134-135] or [5, pp. 67-68, 149]). For ease of computer implementation, the corresponding sampling algorithm is given in pseudocode in Section “Algorithms for sampling”.

Mathematical notation

The most explicit mathematical notation is often also the most cumbersome; this is especially true when dealing with multilevel nesting. We have chosen to use a convention that greatly simplifies the notation of formulas. This convention always and exclusively assigns the same index to a given nesting level. Thus, in this article, the indices h, i, and j correspond to the three levels of nesting, i.e., h for strata, i for PSUs and j for SSUs (Fig. 2). If it is necessary to collapse strata, a fourth level of stratum groups is introduced and associated with the index g (see Section Collapsed strata technique).Fig. 2 Schematic representation of a small example of a hierarchical structure corresponding to populations sampled by stratified two-stage sampling (PSUs in strata and SSUs in PSUs). The indices h, i, and j correspond to the three nested levels, i.e., h for strata, i for PSUs and j for SSUs. A finite population of M PSUs is partitioned into H strata. In turn, a finite population of N SSUs is partitioned by the M PSUs. In the figure, the units are numbered in an ordered fashion to avoid the need for crossing arrows symbolizing the membership relationships.

Fig. 2

A given symbol always and exclusively represents a variable of the same type or with the same meaning, considered at the level associated with a given index. For example, if y denotes the variable of interest, Yg, Yh and Yi correspond to the totals of the variable at the level of group g, stratum h, and PSU i, and Yj corresponds to the value of the variable for SSU j. The same symbol Y (without an index) represents the total at the population level.

We consider finite populations of M PSUs partitioned by strata h=1,2,…,H and of N SSUs partitioned by PSUs i=1,2,…,M (Fig. 2). The uppercase letters M and N represent the numbers of PSUs and SSUs, respectively, at the population and subpopulation levels, whereas their lowercase counterparts refer to the sample sizes. The combination of one of these letters with the index h or i denotes the number of PSUs or SSUs in any of the higher nesting levels. For example, Mh denotes the number of PSUs in stratum h, Ni is the number of SSUs in PSU i, and Nh is the number of SSUs in stratum h (via the PSUs to which they belong).

A sample selected without replacement can be represented in several ways. If the order of the draws does not matter, the sample can be considered an unordered set of labels and represented as an array (or list) of units, e.g., (7,2,3,6,5) for a sample of n=5 out of N=7 units. Owing to the nested structure of the strata, PSUs and SSUs, this choice leads to the simultaneous use of multiple indices; i.e., in this article, h is used for the strata, i for the PSUs and j for the SSUs. An additional index g should be added for stratum groups. Such multiple-index notation is the most common in the literature (e.g., [9, Sec. 8.11] or [10, Sec. 11.7]), but it is not the most convenient from an operational perspective. Another option is to denote the sample as an indicator vector I=(I1,I2,…,Ii…,IN); i.e., for the example above, I=(0,1,1,0,1,1,1). This notation, introduced by Cornfield [11], is very useful for both mathematical derivations and computer implementations. We let Ii (i=1,2,…,M) be the indicator of sample inclusion at the first stage of sampling:(1) Ii={1ifPSUibelongstothefirst-stagesample0otherwise

Ij (j=1,2,…,N) is the second-stage indicator:(2) Ij={1ifSSUjbelongstothesecond-stagesample0otherwise

Regarding the partitioning of a set of J elements into K classes (e.g., partitioning the SSUs by the PSUs or partitioning the PSUs by the strata), at least three simple (equivalent) representations are possible [12, Fig. 1]: (i) a membership vector C indicating the class of each element; (ii) a (J×K) binary matrix Z where each column is a membership indicator corresponding to a class; and (iii) a composite tabular data structure with an indirection vector V for the elements and a (K×2) table A that specifies, for each class, the first and last entries in the vector V. One can use any representation depending on the context. In writing formulas, we use representation (ii). We let Zhi (h=1,2,…,H; i=1,2,…,M) be the membership indicator:(3) Zhi={1ifPSUibelongstostratumh0otherwise

Zij (i=1,2,…,M; j=1,2,…,N) is the following indicator:(4) Zij={1ifSSUjbelongstoPSUi0otherwise

Obviously, Ii=0 implies that Ij=0 when Zij=1. In the design-based framework, the sample indicators Ii and Ij are random variables. Conversely, the indicators Zhi and Zij are fixed because they correspond to the hierarchical structure of the sampled populations (PSUs in strata and SSUs in PSUs).

Parameters

We let the variable of interest be y. The values measured or observed on the SSUs are indexed as Yj (j=1,2,…,N). The totals or means of y are defined by moving up the hierarchy of levels, i.e., through the PSUs and strata, up to the population level:(5) Yi=∑j=1NZijYj=NiY‾i

(6) Yh=∑i=1MZhiYi=NhY‾‾h

(7) Y=∑h=1HYh=NY‾‾

where Y‾i, Y‾‾h and Y‾‾ are the means per SSU at the PSU, stratum and population levels, respectively.

Estimators

For a particular multistage sampling design, such as stratified two-stage sampling with SRSWOR at both stages, one option is to consider explicit formulas, as we do in this article. Another option is to refer to a general implicit formulation, as we did in a previous article [7].

Without considering possible nonresponse (another indicator variable should be introduced to indicate unit nonresponse), we obtain a sample value whenever Ij=1; otherwise, the value remains unknown. This also requires that we have Ii=1 when Zij=1. The formulation of the estimators that we consider in this article involves introducing the sample inclusion indicators Ii and Ij into the summation terms for all units of each population (i.e., PSUs and SSUs). This approach is actually much more convenient than calculating sums directly on samples, as is usually presented in the survey sampling literature.

Estimating a total

The expansion estimator of the total is computed in three bottom-up steps: (i) estimating the total for each PSU i such that Ii=1 (Eq. (8)); (ii) estimating the total for each stratum h=1,2,…,H (Eq. (9)); and finally, (iii) summing the estimates over the strata (Eq. (10)). This leads to Algorithm 3 (Section “Total estimation”).(8) Y^i=Nini∑j=1NIjZijYj=NiY‾^i

(9) Y^h=Mhmh∑i=1MIiZhiY^i=NhY‾‾^h

(10) Y^=∑h=1HY^h=NY‾‾^

The sampling variance of Y^ is estimated in an unbiased way by combining the variance estimators for the SRSWOR of the SSUs from each PSU i such that Ii=1 (Eq. (11)) with those of the SRSWOR of the PSUs from each stratum h=1,2,…,H (Eq. (12)) and finally summing over the strata (Eq. (13)):(11) V^(Y^i)=Ni(Ni−ni)Si2ni

(12) V^(Y^h)=Mh(Mh−mh)Sh2mh︸V1+Mhmh∑i=1MIiZhiV^(Y^i)︸V2

(13) V^(Y^)=∑h=1HV^(Y^h)

where Si2 is the unbiased estimator of the variance for the SSUs in each PSU i such that Ii=1:(14) Si2=1ni−1∑j=1NIjZij(Yj−Y^iNi)2

and where Sh2 is an estimator of the variance for the PSUs in each stratum h=1,2,…,H:(15) Sh2=1mh−1∑i=1MIiZhi(Y^i−Y^hMh)2

This translates into Algorithm 4 (Section “Total estimation”).

Notably, Sh2 is not an unbiased estimator of the variance for the PSUs in stratum h. An unbiased estimator is given as follows:(16) Sh2−1mh∑i=1MIiZhiNi(Ni−ni)Si2ni

In stratum h, the sampling variance of the total estimator Y^h is the sum of two variance components, i.e., the between-PSU variance Vb and the within-PSU variance Vw:(17) V(Y^h)=Vb+Vw

Notably, V1 is not an unbiased estimator of Vb, and V2 is not an unbiased estimator of Vw. Note also that V1 overestimates Vb. For unbiased estimators, see, e.g., [13, p. 181, first remark] for SRSWOR at both stages and [4, p. 137, Result 4.3.1] for the general case, which is not covered in this article.

Estimating a mean

When the parameter of interest is a mean or, equivalently, a proportion when y is a binary variable, we obtain the estimators at different levels from the corresponding estimators for the totals (see Eqs. (8), (9) and (10)). An unbiased estimator of the sampling variance for the overall mean per PSU is as follows:(18) V^(Y‾‾^)=N−2V^(Y^)

For the direct formulation of the sampling variance of Y‾‾^ (i.e., without going through the estimator of the sampling variance of Y^ as in Eq. (18)), we refer the reader to the following.

When the parameter of interest is a mean, we obtain:(19) Y‾^i=Y^iNi=1ni∑j=1NIjZijYj

(20) Y‾‾^h=Y^hNh=MhNhmh∑i=1MIiZhiNiY‾^i

(21) Y‾‾^=Y^N=1N∑h=1HNhY‾‾^h

Formulas (19) to (21) are translated into Algorithm 5 (Section “Mean estimation”). The sampling variance estimator is computed in three steps as follows:(22) V^(Y‾^i)=Ni−niNiSi2ni

(23) V^(Y‾‾^h)=(MhNh)2Mh−mhMhSh2mh+MhNh2mh∑i=1MIiZhiNi2V^(Y‾^i)

(24) V^(Y‾‾^)=1N2∑h=1HNh2V^p(Y‾‾^h)

with(25) Si2=1ni−1∑j=1NIjZij(Yj−Y‾^i)2

and(26) Sh2=1mh−1∑i=1MIiZhi(NiY‾^i−Y‾^h)2

where Y‾^h=NhY‾‾^h/Mh is an estimator of the mean per PSU in stratum h. Formulas (22) to (26) are translated into Algorithm 6 (Section “Mean estimation”).

Estimating a ratio

In addition to a total or a mean, it is often necessary to estimate a ratio R=Y/Z when z is a second variable of interest. For estimating a ratio, a distinction can be made between a separate ratio (or stratum-by-stratum ratio) estimator, which is defined as a weighted sum of stratum ratios, and a combined ratio (or across-stratum ratio) estimator (see, for example, [14, Ch. III, Sec. 6.1], [15, Sec. 5.2] or [16, pp. 46-48]). In this article, we refer to the latter, denoted as R^=Y^/Z^.

Although Y^ and Z^ are both unbiased estimators, unbiasedness is not preserved when forming R^ since a ratio is not a linear combination. Usually, the relative bias is negligible (for details, see [4, Sec. 5.6]). The estimator R^ is consistent in the sense that R^=R when the sample size is equal to the population size (a census), and this property makes sense (finite population consistency; see [3, p. 21]).

Estimating R requires applying Algorithm 3 twice, with y and z as arguments for the variable of interest, to obtain Y^ and Z^. The sampling variance of R^ is approximated by Taylor linearization via the method of Woodruff [1], which greatly simplifies the problem by avoiding the need to consider covariances [6, Sec. 6.5]. For an illustration of this approach in two-stage sampling, see [3, pp. 311-312] or [4, pp. 178-180]; the basic approach can be found in [2, Ch. 7, Sec. 1, 2] or [17, Sec. 7.6], for example. Given Xj=Yj−R^Zj for j such that Ij=1 (i.e., for the sampled SSUs), the sampling variance of R^ is simply estimated by successively applying Algorithms 3 and 4 with x as the variable of interest and finally dividing the overall variance estimate by Z^2.

Despite its practical importance, there is no generally accepted terminology for this method; it is termed, for example, the Tepping-Woodruff method by [4, pp. 497, 499] (after [18] and [1]) and the linear substitute method by [19, pp. 213, 404]. In fact, the Woodruff method is not mentioned as often as one might expect in books on sampling surveys. In particular, this method is used in SAS/STAT software [20] (for an overview of variance estimation methods used in statistical software, we refer interested readers to Carlson [21], for example).

Collapsed strata technique

When using stratified sampling designs, one possible technique is to use many strata and select only one unit per stratum (mh=1 for h=1,2,…,H). This is called one-per-stratum sampling, which, of course, immediately poses a problem for estimating the sampling variance. The same nonestimability problem can occur as a result of nonresponse (unit nonresponse). In this case, it is even possible to obtain mh=0 for some strata. When mh<2 for one or more strata, the sampling variance cannot be estimated via the usual formula (see, for example, Eq. (15)). To overcome this obstacle, one option is to group strata together at the estimation stage to obtain sample sizes that guarantee that mg≥2, where g is the index used for the groups of strata. This technique is called collapsed strata estimation. However, the usual formula for variance estimation cannot be used because grouping strata at the estimation stage does not match the actual structure of the sampling design used to select the sample. Using notation that facilitates computer implementation, we recall the formulas for the estimation when multiple strata are collapsed.

In the literature, the classic situation is when strata for which mh=1 are grouped in pairs (e.g., [10, p. 82]). In what follows, we consider the most general situation, where some strata are not grouped and are treated as usual, while others are grouped to ensure a sample size of at least mmin units (usually, mmin=2 is chosen, but a higher value can also be used). We denote the set of ungrouped strata as H. Therefore, we have mh≥mmin for all h∈H. A set of groups g=1,2,…,G defines a partition of the strata such that h∉H, where Lg≥2 is the number of strata in group g.

Let Zgh (g=1,2,…,G; h=1,2,…,H) be the membership indicator:(27) Zgh={1ifstratumhbelongstogroupg0otherwise

For g=1,2,…,G, we have Zgh=0 for h∈H, and(28) Lg=∑h=1HZgh

Estimation of the total

The population total can be written as follows:(29) Y=∑h=1HYh=∑h∈HYh+∑g=1GYg=∑h∈HYh+∑g=1G∑h=1HZghYh

and is estimated without bias by(30) Y^=∑h=1HY^h=∑h∈HY^h+∑g=1GY^g=∑h∈HY^h+∑g=1G∑h=1HZghY^h

Estimation of the sampling variance

The sampling variance of the estimator (30) can be estimated via a mixed strategy as follows [6, p. 55, Eq. (2.5.9)]:(31) V^mix(Y^)=∑h∈HV^(Y^h)︸VH+∑g=1GLgLg−1∑h=1HZgh(Y^h−AhAgY^g)2︸Vcs

where the Ah values are the values of a size variable defined at the stratum level and where Ag is defined as(32) Ag=∑h=1HZghAh

In the two-term estimator (31), the Ah values are used only for the Vcs term, i.e., for h∉H. However, in what follows, we consider the relationship between Ah and Yh for h=1,2,…,H because there is usually no reason to assume that if the relationship holds for h∉H, it does not hold for h∈H. The expression for Vcs [3, p. 139, Eq. (5A.57)] was introduced by Hansen et al. [2, p. 218, Eq. (5.2)] in the context of estimating a ratio (see also [22, Ch. 9, Sec. 15, 28]). Note that the term collapsed strata technique is not used by Hansen et al. [2], [22]. In addition to being biased, Vcs is inconsistent (in the sense of finite population consistency) [6, p. 53], which means that it has two fundamental drawbacks.

Hansen et al. [2, p. 219] and Wolter [6, p. 51] assumed that Ah values are highly correlated with Yh values, whereas Cochran [3, p. 139] assumed that they are good predictors (which is not the same). Wolter [6, p. 52] noted that, in many real applications, Ah is simply taken to be the number of units within a stratum h. If Ah=1 for h=1,2,…,H, then the expression for Vcs simplifies to [3, p. 139, Eq. (5A.56)], [6, p. 53, Eq. (2.5.8)](33) ∑g=1GLgLg−1∑h=1HZgh(Y^h−Y^gLg)2

The two-term estimator (31) can be used for any stratified design, as long as samples are selected independently in each stratum. In particular, in the case of stratified two-stage random sampling without replacement, V^(Y^h) is given by the expression (12), and Y^h is given by the expression (9). The bias of the two-term estimator comes from the bias of the right-hand term Vcs since the left-hand term VH corresponds to the unbiased estimator for the strata h∈H. Thus, the larger the set H is, the less the two-term estimator (31) is biased.

Formula for the bias of Vcs

To simplify the notation, we use Vh≡V(Y^h). For example, with stratified simple random sampling without replacement, we have:(34) Vh={Mh(Mh−mh)Sh2mhifmh>00otherwise

The sampling variance for Y^g is denoted as(35) Vg=∑h=1HZghVh

The bias of the right-hand term (Vcs) in the two-term variance estimator (31) can be written as follows [6, p. 53, Theorem 2.5.1]:(36) Bias(Vcs)=∑g=1GVAg2−2VAg,VgLg−1Vg︸term1+∑g=1GLgLg−1∑h=1HZgh(Yh−AhAgYg)2︸term2

where(37) VAg2=LgAg2∑h=1HZghAh2−1

and(38) VAg,Vg=LgAgVg∑h=1HZghAhVh−1

In Eq. (36), if Ah=1 for all h, then term 1 vanishes, whereas if Ah=αYh, then term 2 vanishes (see also [6, p. 54]). For the first case, Rust and Kalton [23] performed a thorough investigation of the (relative) bias and concluded that there is no general solution for optimal collapse.

A commonly assumed condition for the bias (or relative bias) of Vcs to be low is that the Ah values are approximately proportional to the Yh values, which implies a high correlation. In this context, we would like to draw the reader’s attention to an aspect of the problem that seems to have been overlooked in the literature. While the correlation is invariant to any linear transformation, this is not true for term 1 in the bias expression (Eq. (36)). As a result, the magnitude of the bias can be very different for two situations that are otherwise identical in terms of the correlation between the Ah values and the Yh values. In our view, this is another drawback for Vcs because even if one can guarantee a high correlation between the totals of the size variable and those of the variable of interest at the stratum level, this does not ensure that the bias (or relative bias) of Vcs will be low. Although, on average, the bias decreases with increasing correlation, a high correlation does not ensure low bias. Importantly, the Ah values should be good surrogates for the Yh values — not just highly correlated — which seems very difficult to guarantee in practice. Therefore, in general, we believe that using Vcs for variance estimation with collapsed strata involves taking a gamble on the bias, which is a situation that we find unsatisfactory. This does not mean that the bias cannot be small, only that we have no way of ensuring that this is indeed the case. It is beyond the scope of this article to explore this issue further or to consider alternative solutions, which may be the subject of another article.

Algorithms

The proposed algorithms (procedures and functions) are presented as pseudocodes with classical control structures. There are three possible passing modes for a procedure parameter depending on whether it is only for input (”in” mode, where its value remains unchanged), only for output (”out” mode, where its value is defined within the procedure) or for both input and output (”in–out” mode, where its value is modified within the procedure). We use a down arrow, an up arrow or a down–up arrow above the symbol of a parameter to indicate the ”in”, ”out”, and ”in–out” passing modes, respectively. By definition, all the parameters of a function are in input-only passing mode; therefore, the downward arrows are omitted.

Since stratified two-stage sampling involves three levels — it is a particular instance of three-stage sampling — the total numbers of strata, PSUs and SSUs are denoted as N1, N2 and N3, respectively. When collapsing strata is necessary, the number of groups is denoted as N0 for convenience. The values of the variable of interest are in an (N3×1) array Y3, the total or mean estimates are in an (N2×1) array Y2 at the PSU level and in an (N1×1) array Y1 at the stratum level, and the estimate Y is at the population level.

In the algorithms, the partitions are represented by one-dimensional membership arrays replacing the binary matrices (see Section Mathematical notation), namely, an (N3×1) array C3 for the partitioning of the SSUs by the PSUs and an (N2×1) array C2 for the stratification of the PSUs. Similarly, when collapsing strata is necessary, the partitioning of the strata by groups is represented by an (N1×1) array C1, where ungrouped strata have a zero group code. The variance estimates v1 and v2 are intermediate quantities that do not need to be stored in an array; therefore, they appear without indices in our algorithms. To limit the number of arguments in Algorithms 2, Algorithm 3, Algorithm 4, Algorithm 5 and 6, the cardinalities N2 and N3 are derived from the associated membership vectors C2 and C3 using the Card(·) operator (this is not mandatory, of course). When they are not identical, Table 1 presents the correspondences between the algorithmic and mathematical notations.Table 1 Table of correspondences between algorithmic and mathematical notations.

Table 1	Algorithmic	Mathematical	
Number of groups (of strata)	N0	G	
Number of strata (of PSUs)	N1	H	
Number of PSUs	N2	M	
Number of SSUs	N3	N	
Number of SSUs in stratum h	NHh	Nh	
Estimate at the population level	Y	Y^ or Y‾‾^	
Estimate at the h-th stratum level	Y1h	Y^h or Y‾‾^h	
Estimate at the i-th PSU level	Y2i	Y^i or Y‾^i	
Value of the variable of interest for SSU j	Y3j	Yj	
Group membership for stratum h	C1h	Zgh	
Stratum membership for PSU i	C2i	Zhi	
PSU membership for SSU j	C3j	Zij	
First-stage sample inclusion indicator for PSU i	I2i	Ii	
Second-stage sample inclusion indicator for SSU j	I3j	Ij	
Variance estimate for the PSUs in stratum h	v1	Sh2	
Variance estimate for the SSUs in PSU i	v2	Si2	
Sampling variance estimate at the population level	V	V^(Y^) or V^(Y‾‾^)	
Sampling variance estimate at the h-th stratum level	V1h	V^(Y^h) or V^(Y‾‾^h)	
Sampling variance estimate at the i-th PSU level	V2i	V^(Y^i) or V^(Y‾^i)	

Algorithms for sampling

In this section, we consider situations where the populations can be stored in random access memory (RAM). Therefore, we do not consider sampling files with sequential access or other types of data streams.

Stratified two-stage sampling requires a procedure for selecting units within groups that form a partition to sample (i) PSUs in strata (first stage of sampling) and (ii) SSUs in PSUs selected at the first stage (second stage of sampling). Formally, the two stages are identical, except that the strata are not sampled, unlike the PSUs. Therefore, the same procedure can be used, with a logical array indicating whether a group (stratum and then PSU) is involved in the subsampling process. With a suitable data structure, it is then sufficient to make two successive calls to this procedure to obtain a stratified two-stage sample as the output. The same approach can be used for more levels, but it is also possible to rely on the data structures and algorithms described in [7].

Sampling within groups of units that form a partition

We consider a population of units of indices j=1,2,…,N and groups i=1,2,…,M that form a partition of such units (e.g., SSUs in PSUs). To browse the units of a group, we refer to a composite data structure with an indirection (N×1) array B for the units and an (M×2) array A that contains, for each group, the first and last entries in array B [12, Fig. 1].

We define an (M×1) array I1 such that I1i is TRUE if group i is involved in the subsampling process and FALSE otherwise. Here, we consider that for each group i for which I1i is TRUE, ni out of Ni units have to be selected by SRSWOR. We denote the algorithm for sampling n0 out of N0 units with equal probability and without replacement asSRSWOR(N0↓,n0↓,I0↑)

This algorithm returns a sample membership indicator I0 such that I0k is TRUE if unit k belongs to the sample and FALSE otherwise. The choice of the SRSWOR algorithm is left to the discretion of the reader (e.g., see [5, Sec. 3.9]).

Algorithm 1 selects samples for all groups i such that I1i is TRUE and returns an (N×1) array I2 such that I2j is TRUE if unit j belongs to the sample and FALSE otherwise.Algorithm 1 Sampling within groups of units that form a partition.

Algorithm 1

Stratified two-stage simple random sampling without replacement

In this section, we assume that the programming language that will be used for computer implementation has a mechanism for preserving the values of some variables between two calls of a procedure, i.e., either via a static declaration of local variables (e.g., in C/C++) or by using global variables (e.g., in Pascal/Delphi). This is not mandatory, but the idea is to encapsulate the preprocessing of the composite data structure inside the procedure, thereby simplifying the use of the procedure and reducing the number of arguments, while avoiding the need to repeat this initialization when multiple calls are needed, e.g., in the case of a Monte Carlo simulation. Thus, a logical variable F (i.e., a flag) is set to TRUE before the first call of the procedure described by Algorithm 2. This procedure can be easily extended to more sampling levels. We refer the reader to [12, p. 5] for the definition of the following procedure:BuildTableOfClasses(C↓,N↓,H↓,A↑,B↑)

Algorithm 2 Stratified two-stage sampling with SRSWOR at both stages.

Algorithm 2

Algorithms for estimation

Total estimation

Formulas (8) to (10) for estimating the totals are translated into Algorithm 3. Formulas (11) to (15) for estimating the corresponding sampling variances are translated into Algorithm 4, which assumes that mh≥2 for h=1,2,…,N1 and that ni≥2 for i=1,2,…,N2.Algorithm 3 Estimating the totals for the sampled PSUs (Y2), strata (Y1), and population (Y).

Algorithm 3

Algorithm 4 Estimating the sampling variances for the totals in Y2, Y1, and population total Y.

Algorithm 4

Mean estimation

Formulas (19) to (21) for estimating the means are translated into Algorithm 5. Formulas (22) to (26) for estimating the corresponding sampling variances are translated into Algorithm 6, which assumes that mh≥2 for h=1,2,…,N1 and that ni≥2 for i=1,2,…,N2.Algorithm 5 Estimating the means for the sampled PSUs (Y2), strata (Y1), and population (Y).

Algorithm 5

Algorithm 6 Estimating the sampling variances for the means in Y2, Y1, and population mean Y.

Algorithm 6

Algorithms for collapsing strata

In the case where mh<2 for one or more strata, the sampling variance cannot be estimated with the usual estimator. The default option is to assign a zero value to the variance estimate in these strata. SAS/STAT software does this when the NOCOLLAPSE option of the SURVEYREG procedure is used [20, p. 10134]. Another option is to collapse strata. This operation is often performed manually, but we want to be able to perform it automatically, particularly so that it can be examined in Monte Carlo simulations. Therefore, we must determine a set of rules that can be translated into a procedure that is guaranteed to end in a reasonable computing time. For example, in the case of the SURVEYREG procedure in SAS/STAT software, the rules are very simple [20, p. 10134]:• If there are multiple strata that each contain only one sampling unit, then all these strata are combined into a single group.

• If there is only one stratum containing only one sampling unit, then this stratum is aggregated with the previous stratum (the strata are ordered). If this is the first stratum, then it is aggregated with the next stratum.

Other rules can be used (see, for example, the function collapse.strata in the R package ReGenesees, [24]), but they must allow computation to be performed without human intervention.

In this section, we present two algorithms for collapsing strata — one deterministic and the other stochastic — depending on the information available for grouping the strata. Note that this must not be information from the sample at hand, which would lead to a severely downward biased estimator of the sampling variance (see [3, p. 139], [23, p. 71], [6, p. 54] or [10, p. 82]). Specifically, our two algorithms are based on how the Ah values are defined.

The inputs to both algorithms are (1) the number of strata N1; (2) the initial sample sizes in an (N1×1) array m; and (3) the limit value min for the number of units in a stratum or group of strata. The outputs are (4) the number of groups N0; (5) an (N1×1) array C1 indicating the membership of the strata in the groups (C1 is the algorithmic counterpart of Zgh); and (6) an (N1×1) array L containing the sizes of the groups (the number of groups N0 is known only afterwards; therefore, the default size of this array is N1).

The principle of our collapsing strata algorithms is very similar to that of agglomerative hierarchical clustering, except that the binary tree of agglomerative nodes is fully constructed only if all strata are to be aggregated together. In the general case, the result of the algorithm is a partition of the strata into several groups, with each group corresponding to a subtree. Strata that are not aggregated with others remain singletons, and their group code is zero (these strata we denote as h∈H).

Our algorithms use a two-column array D to indicate the left (first column) and right (second column) descendants for each node resulting from the agglomeration process. The current size of D is denoted as d. In the case where all strata must be aggregated, the algorithms terminate with d=2×N1−1 since there are N1 leaves and N1−1 nodes. Our algorithms also use (i) an array I to recycle the entries from 1 to N1 to be used for the groups to be formed and (ii) an array S to store the number of units.

Algorithms 7 and 8 are self-explanatory and need no further introduction. Algorithm 9 describes the final step in building groups from agglomerative nodes. It consists of identifying the nodes that are the starting points for forming the groups and deploying them by simulating recursion via a stack. Internally, it uses an indicator array I to identify the starting-point nodes to be deployed to form groups.Algorithm 7 Initializing data structures for aggregation.

Algorithm 7

Algorithm 8 Aggregating sets of strata a and b by creating a new agglomerative node.

Algorithm 8

Algorithm 9 Forming groups of strata from subtrees.

Algorithm 9

Cases where strata have unit Ah values for grouping

If Ah=1 for h=1,2,…,H, we assume that the strata are ordered; that is, the numerical indices of the strata define an ordinal variable with H categories such that Y‾1<Y‾2<⋯<Y‾h<⋯<Y‾H. Therefore, two consecutive strata are assumed to be more similar to each other than to other strata. We use the term adjacent strata to refer to two strata h and h′ such that |h−h′|=1. In this case, our proposed algorithm consists of grouping adjacent strata, and if there are two candidate strata, the stratum containing the fewest sampled units is selected. This process is repeated until all stratum groups have a sample size of mg≥mmin. Grouping multiple strata in an agglomerative manner involves updating adjacency relationships, which is reflected in Algorithm 10. Importantly, the update performed is correct only if the strata are aggregated sequentially. A change in the way strata and stratum groups are aggregated must involve a change in Algorithm 10. The proposed stratum aggregation procedure is deterministic and translates into Algorithm 11. The general structure of this algorithm is shown in Fig. 3.Algorithm 10 Updating adjacency relationships between sets of strata after a and b were aggregated.

Algorithm 10

Algorithm 11 Strata aggregation based on the relationship of adjacency.

Algorithm 11

Fig. 3 Flowchart of Algorithm 11. The term set of strata designates either a single stratum (singleton) or a group of strata.

Fig. 3

Cases where strata have different Ah values for grouping

When Ah is a size value, we first consider the special case where Ah=αYh, that is, when the Ah values are strictly proportional to the Yh values (h=1,2,…,H); this implies a perfect correlation of the size variable with the variable of interest at the stratum level. In this case, the best stratum grouping will minimize the expression for the bias (Eq. (36)) to (almost) zero. Recall that when Ah=αYh, term 2 in the bias expression is zero, which means that the active term in the minimization is term 1, which translates into Algorithm 12.Algorithm 12 Term 1 in the bias of Vcs.

Algorithm 12

In the general case, the Ah values are at best only approximately proportional to the Yh values, which is formally written as Ah=αhYh, with a variable coefficient αh for h=1,2,…,H; this results in a variable correlation between the Ah values and the Yh values. In this case, the best stratum grouping minimizes the heterogeneity of the Ah values within each group of strata. In this article, heterogeneity is defined in the sense of inertia. The total inertia T is partitioned as follows:(39) T=∑h=1H(Ah−A‾)2=∑g=1G∑h=1HZgh(Ah−A‾)2=∑g=1G∑h=1HZgh(Ah−A‾g)2︸W+∑g=1GLg(A‾g−A‾)2︸B

with the overall mean(40) A‾=1H∑h=1HAh

and the intragroup mean (for g=1,…,G)(41) A‾g=AgLg=1Lg∑h=1HZghAh

and where W and B are the within- and between-group inertia, respectively. Hence, in this article, minimizing the heterogeneity of the Ah values within each group of strata involves minimizing W, or equivalently, maximizing B (since T is constant). The objective function W to be minimized is translated into Algorithm 13.Algorithm 13 Within-group inertia.

Algorithm 13

Obtaining stratum groups by minimizing a given objective function involves solving a combinatorial optimization problem. In what follows, we seek to obtain only a good solution, not necessarily the global optimum. To achieve this goal, the simplest heuristic is to explore the space of possible solutions via local (stochastic) improvement of initial completely random solutions. The search process repeats the following sequence of three operations several times: (1) generating an initial solution by random grouping; (2) evaluating the objective function for the initial solution; and (3) improving the initial solution. The best grouping stored during this process is returned as the output.

The procedure corresponding to operation (1) is to randomly select one of the strata containing too few sampled units and merge it with another randomly selected stratum. The process is repeated as long as there are strata (or groups of strata) containing too few sampled units. The random grouping of strata is translated into Algorithm 14. This procedure outputs groups of strata that satisfy the sample size constraint specified by the value of min and is independent of the objective function. The evaluation of the objective function for the initial random solution in operation (2) ends the initialization for the current iteration of the search process. The local stochastic improvement of the initial solution in operation (3) consists of randomly modifying the initial solution if possible while maintaining the constraint on the number of units in each group and avoiding group degeneracy (a group necessarily contains at least two strata). This improvement procedure stops when a certain maximum number of iterations is reached with no improvement; this indicates that a local minimum of the objective function has been reached. In the case of an apparent deadlock (there is no stratum that can change groups, or there is no group that can receive a stratum), the procedure simply stops and proceeds to another initial solution. This procedure is described in Algorithm 15, which uses Algorithm 16, whose role is simply to update the group sizes and the numbers of units they contain.Algorithm 14 Grouping strata in a random way.

Algorithm 14

Algorithm 15 Stochastic local improvement of an initial grouping of strata.

Algorithm 15

Algorithm 16 Updating the size and number of units in a (which receives) and b (which gives).

Algorithm 16

The resulting heuristic involving operations (1), (2) and (3) translates into Algorithm 17, where max is the maximum finite value for the floating-point type used in the computer implementation of the algorithm (i.e., with single, double, extended or quadruple precision). The general structure of this algorithm is shown in Fig. 4. Algorithm 18 is self-explanatory and needs no further introduction.Algorithm 17 Strata aggregation based on a size variable defined at the stratum level.

Algorithm 17

Fig. 4 Flowchart of Algorithm 17. The term set of strata designates either a single stratum (singleton) or a group of strata. The number of allowed iterations is K.

Fig. 4

Algorithm 18 Saving or restoring a grouping and its associated objective function value.

Algorithm 18

In practice, the goal is to minimize the within-group inertia (Algorithm 13). As an objective function, term 1 of the bias can only be used in the case of a known population, which also implies that Vh is known for h=1,2,…,H. These conditions make it possible to check that the heuristic minimizes the bias to (almost) zero when Ah=αYh for h=1,2,…,H. We use this approach in Section “Method validation” with the following modifications to Algorithm 17: (i) in line 1, V1 is added to the parameter list of the procedure; (ii) in line 8, Algorithm 12 is used as the objective function instead of Algorithm 13; and (iii) in line 10, the comparison is made between absolute values since the bias can be negative [6, p. 54], even though it is usually positive [2, Ch. 9, Sec. 5]. The same changes are applied to Algorithm 15 (lines 1, 44 and 45, respectively).

Regardless of the objective function, in a simulation context, we can compute the bias expression (36) as described in Algorithm 19.Algorithm 19 Bias of Vcs.

Algorithm 19

Method validation

In this section, validation consists of checking that the computer implementations of the algorithms return the correct results. Validating the algorithms themselves would require formal proofs, a subject far beyond the scope of this article.

Data example

The illustrative example in this section may be useful for readers in checking their own computer implementations in the programming language of their choice. The example dataset is provided in Section “Supplementary material”.

We consider a simple situation where a spatial domain of interest is bounded by an iso-oriented rectangle (i.e., its sides are parallel to the abscissa and ordinate axes), lying between xmin and xmax on the abscissa and between ymin and ymax on the ordinate. This domain is discretized by a population of N=2500 square cells of equal size, organized according to a grid of NX=100 columns and NY=25 rows. For simplicity, we set xmin=ymin=0, xmax=100 and ymax=25 (Fig. 5.a). The population of the grid cells corresponds to the population of the SSUs. In the dataset, the SSUs are numbered from 1 to N by traversing the grid from top to bottom and from left to right. The values Yj for j=1,2,…,N are obtained by counting the objects in each grid cell. A total of 10000 objects are distributed among the N=2500 grid cells via allocation with unequal probabilities, as explained in [25, Sec. 5.2]. The example used here corresponds to that shown in [25, Fig. 12.1a], with a pronounced gradient along the x-axis, from xmin to xmax (Fig. 5.a). Two PSUs are defined in each column, with the size of the top PSU randomly drawn uniformly between 8 and 15 SSUs (Fig. 5.b and 5.c). Thus, the population of PSUs contains M=200 units of variable size (in terms of the number of SSUs). The strata are defined by partitioning the columns into H classes (e.g., H=5 in Fig. 5.a). This partitioning is performed by minimizing the within-group inertia, which is computed on the column totals (or equivalently, the column means). The exact solution is obtained via the algorithm of Fisher [26], which was efficiently implemented by Hartigan [27, Ch. 6].Fig. 5 Schematic representation of the example dataset. (a) A grid of NX=100 columns and NY=25 rows defines a population of N=2500 secondary sampling units (SSUs). The marginal graph (in gray) represents the mean of the grid cells along the x-axis. This reflects the gradient from 0 to 100. The 100 columns are partitioned in the direction of the gradient. In this example, we have strata h=1,2,…,5 (H=5). (b) Details of the left part of the grid, from columns 1 to 10. (c) Details of the right part of the grid, from columns 90 to 100. In each column, the two PSUs are shown in yellow and green. The PSUs consist of 8 to 17 SSUs.

Fig. 5

Stratified two-stage sampling and estimation procedures

If the sampling procedure and the procedures corresponding to the estimators are consistent, then a Monte Carlo simulation based on them should produce results that match the expected theoretical results. This is a simple, universal approach to checking procedures in the field of probability sampling.

To check that the procedure implementing the sampling design is correct, as a first step, we can compare the (first-order) inclusion probabilities for the PSUs (πi for i=1,2,…,M) and the overall inclusion probabilities for the SSUs (πj for i=1,2,…,N) with the values approximated by the Monte Carlo simulation [28, Sec. 2.3]. The theoretical probabilities are πi=mh/Mh when Zhi=1 (h=1,2,…,H; i=1,2,…,M) and πj=πini/Ni when Zij=1 (i=1,2,…,M; j=1,2,…,N). Their Monte Carlo counterparts are computed from the indicators I2 and I3 (respectively) by repeating the sampling process many times.

If the examination of the inclusion probabilities does not reveal any discrepancies that might call into question the sample selection procedure, the second step is to check the computer implementation of the estimators; this can be done with a Monte Carlo simulation based on many sample selections. For each sample, we can evaluate the estimator of the total (or mean) as well as that of its sampling variance. We know the values of the total (here, 10000 objects) or the mean (here, 4 objects per SSU), which we can compare with the Monte Carlo approximations of the expectations of the corresponding estimators. The estimators of the total or the mean are design unbiased, so we should observe only a very small discrepancy between the parameter values and the Monte Carlo approximations of the expectations of their estimators (depending on the simulation effort). The same applies to the estimation of the sampling variances for the total and mean estimators; the values approximated by Monte Carlo simulation should match those of the approximations of the expectations of the sampling variance estimators. The sampled populations of PSUs and SSUs are fully known in the simulation, so we can also make comparisons with the theoretical sampling variances. Since the theoretical variances cannot be computed in practice, we do not present their formulas or the corresponding procedures in this article. However, it is straightforward to obtain the theoretical sampling variance formulas from those of their estimators, and the same is true for the associated procedures. Notably, the totals or means at the PSU, strata and population levels can also be computed simply by changing the arguments passed to the procedures that compute the corresponding estimates. First, all values of I2 and I3 are set to TRUE. Next, for the total, the procedure call is as follows:

EstimateTotals(N1,C2,I2,M,M,C3,I3,N,N,Y3,Y2,Y1,Y)

Similarly, for the mean, the procedure call is as follows:

EstimateMeans(N1,NH,C2,I2,M,M,C3,I3,N,N,Y3,Y2,Y1,Y)

We use the dataset example with H=5 strata and Mh=104,36,28,20,12 for h=1,2,3,4,5 (Fig. 3.a). We set the sampling effort for both sampling stages and replicate the sample selection 107 times.

For example, to check the inclusion probabilities, in the first stage of sampling, we use a uniform allocation rule with mh=4 for h=1,2,3,4,5, and in the second stage of sampling, we compute ni=⌈Ni/3⌉ for i=1,2,…,M. For both sampling stages, the agreement between the theoretical and approximated inclusion probabilities is indisputable (Fig. 6).Fig. 6 Relationship between the theoretical values and Monte Carlo approximations of the inclusion probabilities. (a) Relationship in the first stage of sampling between the theoretical values (πi) and their Monte Carlo approximations π^i. (b) Relationship in the second stage of sampling between the theoretical values (πj) and their Monte Carlo approximations π^j.

Fig. 6

For example, to check the estimators, in the first stage of sampling, we use the allocation mh=6,4,4,4,2 for h=1,2,3,4,5, and in the second stage of sampling, we use a uniform allocation rule with ni=4 for i=1,2,…,M.

For the estimation of the population total, the procedures to be validated are EstimateTotals and EstimateTotalVariances. The computed values are rounded to the nearest unit. The Monte Carlo approximation for the expectation of the population total estimator is 10000, which indicates that the estimator is design unbiased. The Monte Carlo approximation for the expectation of the sampling variance estimator is 626547, which should be compared with the Monte Carlo approximation of the sampling variance (626750) and the theoretical sampling variance (626488).

Similarly, for the estimation of the population mean, the procedures to be validated are EstimateMeans and EstimateMeanVariances. The computed values are rounded to 6 decimal places. The Monte Carlo approximation for the expectation of the population mean estimator is 4.000082, which indicates that the estimator is design unbiased. The Monte Carlo approximation for the expectation of the sampling variance estimator is 0.100248, which should be compared with the Monte Carlo approximation of the sampling variance (0.100280) and the theoretical sampling variance (0.100238).

The results for the estimation of the total and the mean are as expected, and the two sets of results are also in perfect agreement with each other (the total is divided by N to obtain the mean, and the same is done with N2 for the sampling variance). Thus, these results do not reveal any anomalies in the computations performed. Although this is not a proof, it is reasonable to conclude that the corresponding procedures can be used with confidence.

Stratum aggregation procedures

In any stratum aggregation procedure, the sample size constraint must be satisfied; this immediately leads to the need for a Boolean-valued function that confirms that this is indeed the case (Algorithm 20).Algorithm 20 Checking that the system of strata and groups of strata has valid sample sizes.

Algorithm 20

Deterministic aggregation

For our deterministic stratum aggregation procedure (AggregateStrataDeterministic, Algorithm 11), a specific Boolean-valued function checks whether the adjacency constraint is satisfied (Algorithm 21). This function is correct only if the strata are aggregated sequentially. Algorithm 22 determines whether a list of strata is valid. If a list contains only one stratum, then its group code must be zero. Examples of valid and invalid systems of strata and groups of strata are presented in Table 2.Algorithm 21 Checking that groups are made up of adjacent strata.

Algorithm 21

Algorithm 22 Checking that a list of strata is valid.

Algorithm 22

Table 2 Examples of valid and invalid systems of strata and groups of strata for N1=6.

Table 2	Valid	Invalid	Valid	Invalid	Invalid	Valid	Invalid	
1	0	0	1	0	0	1	0	
2	1	1	1	1	1	1	2	
3	1	1	1	1	1	2	2	
4	0	2	2	0	0	2	0	
5	2	2	2	1	2	3	1	
6	2	3	0	1	3	3	1	

To systematically check procedure AggregateStrataDeterministic, one can generate all sample size allocations for a given number of strata and an overall sample size. We let the overall sample size be n and the number of strata be m. In the context of collapsed strata, n can be decomposed into a sum of m integers such that the order is important and parts with zeros are allowed. The problem to solve is to generate all integer compositions of n into exactly m parts such that(42) (n,0,…,0︸mparts),(n−1,1,…,0︸mparts),…,(0,…,1,n−1︸mparts),(0,…,0,n︸mparts)

If parts with zeros are allowed then the number of integer compositions of n into exactly m parts can be computed as (n+m−1n) [29, Ch. 5, Eq. (2)]. An algorithm that generates integer compositions one by one in the order given in expression (42) can be found in [29, Ch. 5] and is translated into Algorithm 23, where t and h are static local variables (or global variables) such that their values are preserved between two calls. The parameter F is a flag that is set to TRUE before the first call and that eventually takes the value TRUE when the last integer composition has been generated.Algorithm 23 Generating the next integer composition.

Algorithm 23

Finally, the checking procedure for AggregateStrataDeterministic leads to Algorithm 24, where n2 is the overall sample size (e.g., the number of PSUs sampled in stratified two-stage sampling).Algorithm 24 Checking that the deterministic aggregation of strata returns valid outputs for all integer compositions.

Algorithm 24

Stochastic aggregation

For our stochastic stratum aggregation procedure (AggregateStrataStochastic, Algorithm 17), we first check whether the bias is minimized to (almost) zero when the size values are strictly proportional to the total values at the stratum levels (Ah=αYh for h=1,2,…,H, where α is constant). We call this the ideal case.

Next, if the proportionality between the size values and the total values at the stratum levels is only approximate (Ah=αhYh for h=1,2,…,H, where αh is variable), we check whether, on average, the bias decreases as the correlation between the Ah values and the Yh values increases. We call this the approximate case.

In a simulation context, the Ah values can be generated according to a given correlation with the Yh values by drawing random values for αh (see Appendix, Algorithm 25). This approach corresponds to the theme of the variance estimator in the collapsed strata technique since the bias of the estimator Vcs is related to the coefficients αh.Algorithm 25 Using a random multiplication factor to generate a correlated size variable.

Algorithm 25

We use the example dataset with H=10 strata defined by minimizing the within-group inertia (Fisher algorithm), with Mh=60,28,30,22,18,10,12,8,6,6 for h=1,2,…,10. To perform the computations, we use the total at the PSU level as the variable of interest (there is no need to consider the SSUs here). We set the sampling effort in the strata with a proportional allocation of 30 units, that is, mh=9,4,5,3,3,1,2,1,1,1 for h=1,2,…,10.

For the ideal case, the goal is to minimize the bias by using the function Term1Bias (Algorithm 12) as the objective function. We use a constant coefficient α=1, which implies a perfect correlation between the Ah values and the Yh values (the value of α does not matter, as long as it is constant). We call the procedure AggregateStrataStochastic with min=5, K1=K2=100 and the input data shown in Fig. 7a. As the output results, the procedure returns N0=3 stratum groups, as shown in Fig. 7b and 7 c. The relative bias of Vcs with respect to the theoretical sampling variance for the population mean (or total) estimator is 0.000958%, which is almost zero, as expected.Fig. 7 Ideal case (α is constant). The input data and output results are shown for the stochastic aggregation procedure. (a) Input data for h=1,2,…,10. (b) Group codes for h=1,2,…,10. (c) Group sizes for g=1,2,3.

Fig. 7

For the approximate case, the goal is to minimize the within-group inertia by using the function WithinGroupInertia (Algorithm 13) as the objective function. We use randomly generated coefficients αh for the theoretical correlations, which are randomly drawn uniformly between 0.8 and 1. Owing to the small number of strata, the actual correlation between the Ah values and the Yh values fluctuates and is often smaller than the prescribed theoretical correlation. However, below, we present only the relative bias obtained for situations where the actual correlation between the Ah values and the Yh values is greater than 0.8. We call procedure AggregateStrataStochastic with min=5 and K1=K2=100 for 104 sets of Ah values computed by using the αh generated as described above. We compute the relative bias of Vcs with respect to the theoretical sampling variance for the population mean (or total) estimator. We plot the relative bias against the actual correlation between the Ah values and the Yh values (Fig. 8). As expected, the relative bias tends to decrease as the correlation increases. However, for the same correlation value, there can be considerable variability in the corresponding relative bias, and this variability increases as the correlation decreases. This result indicates that the bias of the sampling variance estimator for the collapsed strata technique is not invariant to linear transformations.Fig. 8 Approximate case (αh is variable) using 104 sets of Ah values. Plot of the relative bias (%) against the correlation between the Ah values and the Yh values, for correlations greater than 0.8

Fig. 8

Additional information

Reasons for using a stratified two-stage sampling design

There are at least three practical reasons for using a stratified two-stage sampling design. The first is that it combines the advantages of both stratification and multistage sampling in terms of efficiency and cost. A second reason, specific to stratified sampling, is that it enables estimates to be computed directly at the stratum level without the complications of estimating parameters for unplanned subdomains. A third reason, specific to two-stage sampling, is the possibility of defining the frame of SSUs only for the PSUs selected in the first stage, which can yield considerable economy in terms of preparatory work.

Simplified variance estimators

Simplified variance estimators are not discussed in this article, as we believe that the search for simplification is not relevant now that all computations are performed on (highly efficient) computers. What matters is that the estimators are programmed correctly, regardless of their level of complexity (see, for example, [7]). For an overview of simplified variance estimators, we refer the interested reader to [4, Sec. 4.6], for example.

(Not too) deep stratification

In deep stratification, one or more auxiliary variables correlated with the variable of interest are used to define many strata. In this case, the allocation of overall sampling effort among the strata may result in undersized (or even empty) within-strata samples. At the design stage, one way to address this problem is to use Algorithm 17 to automatically group strata to achieve a minimum sample size in each final group. Minimizing the within-group inertia is an appropriate approach for maximizing the precision of the expansion estimator for a linear parameter (population total, mean or proportion) in the case of stratified simple random sampling. In a univariate context, we simply use Algorithm 13 as is. Adapting this algorithm to a multivariate context is straightforward and is left to the reader.

Algorithm 17 is a convenient way to redefine strata and adapt to any overall sampling effort and allocation rule among strata. In any case, this approach is more practical than iteratively redefining strata from scratch through trial and error until the sample size constraint is satisfied. Another option, which is beyond the scope of the present article, is to rely directly on a tree-based approach for strata formation, as proposed by Benedetti et al. [30], for example. In both cases, the goal is to use automated procedures to reduce the human workload and provide reproductible results. It would be interesting to compare the pros and cons of the two types of approaches, depending on the objectives and data available.

Limitations

None.

Ethics statements

None.

CRediT authorship contribution statement

Philippe Aubry: Conceptualization, Methodology, Software, Writing – original draft.

Declaration of competing interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A

Notation for random variables

Random variables are usually written in upper case, whereas their realizations are written in lower case. Another possibility is to use Dutch notation, which consists of underlining symbols for the random counterparts of fixed quantities [31]. Although this convention for distinguishing between fixed and random quantities is more generally applicable than the previous convention is, it has the disadvantage of making the notation considerably more cumbersome. In this appendix, we choose not to distinguish random variables from their realizations in the notation but simply refer to the context.

Generating a correlated size variable

We let y be the variable of interest and U be a finite population of size N. The set of y values for U is a fixed vector y=(y1,y2,…,yN). We let z be a size variable that is also defined on U, resulting in a vector z=(z1,z2,…,zN) correlated to y with a (finite population) correlation Ryz. One can attempt to directly generate the vector z with the desired correlation Ryz, or one can use a statistical superpopulation model for z and consider a theoretical correlation ρyz. With the second option, since z is a sample of values drawn from a superpopulation, Ryz does not necessarily match ρyz (unless N is very large). In this case, it is possible to generate many instances of the vector z and keep one that gives a value of Ryz that is close to ρyz with respect to a certain level of discrepancy (e.g., when |Ryz−ρyz|<10−6). This is the approach taken in [25, Sec. 4.3]. In what follows, we consider a superpopulation model with a prescribed correlation 0<ρyz<1. In this context, a vector z can be generated from a given vector y and a correlation ρyz defined at the superpopulation level by attenuating the perfect positive correlation. This can be performed by adding a Gaussian error term ϵ∼N(0,σϵ2) (for example) to the y values. The solution for determining the appropriate value of σϵ for obtaining a given value of ρyz can be found in [25, Eq. (9)]. Then, z can be mapped to any given variation range since the correlation is invariant to linear transformations.

Instead of adding an error term to the y values, another approach is to multiply the y values by a random coefficient α that varies in an appropriate way. Thus, the size variable is generated as zi=αiyi for i∈U. For the extreme case where ρyz=1, we have αi=1 for all i∈U. To obtain a correlation ρyz<1, we simply vary the α values so that 1−ϵ<α<1+ϵ with ϵ>0, according to a uniform probability distribution (in this case, we have E(α)=1). Therefore, the larger ϵ is, the farther from 1 the correlation ρyz is. The question then becomes how to determine the appropriate value of ϵ for obtaining a given value of ρyz. From [32, Eq. (12)], we have the following general result:(43) Cov(y,αy)=E(α)V(y)+E(y)Cov(α,y)+E{Δ(α)[Δ(y)]2}

with Δ(α)=α−E(α) and Δ(y)=y−E(y). The two random variables α and y being stochastically independent implies that they are uncorrelated, i.e., the central cross moments are zero; here, Cov(α,y)=0 and E{Δ(α)[Δ(y)]2}=0. Therefore, expression (43) simplifies to Cov(y,αy)=E(α)V(y). Moreover, since we have E(α)=1, we obtain Cov(y,αy)=V(y). Furthermore, since α and y are stochastically independent, we have the following [33, Eq. (2)] [32, Eq. (9)]:(44) V(αy)=V(α)V(y)+V(α)[E(y)]2+V(y)[E(α)]2

Therefore, with z=αy and E(α)=1, the correlation of interest is as follows:(45) ρyz=Cov(y,αy)V(y)V(αy)=V(y)V(y){V(α)V(y)+V(α)[E(y)]2+V(y)}

With ρyz>0 we can write:(46) V(α)=1ρyz2−11+[E(y)]2V(y)

As expected, ρyz=1 implies V(α)=0 and therefore αi=1 for i∈U or, equivalently, zi=yi for i∈U. With α uniformly distributed over the interval ]1−ϵ,1+ϵ[, we have V(α)=ϵ2/3. Thus, after calculating V(α) (Eq. (46)), we obtain the value ϵ=3V(α), which corresponds to the correlation ρyz. This translates directly into Algorithm 25, where MY and VY are the mean and variance, respectively, for the values in the vector Y; R is the correlation ρyz; and Uniform is a function that returns a pseudorandom number on the interval ]0,1[.

Supplementary material

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.mex.2024.102928

Appendix B Supplementary materials

Supplementary Data S1

Supplementary Raw Research Data. This is open data under the CC BY license http://creativecommons.org/licenses/by/4.0/

Supplementary Data S1

Data availability

I have shared the example dataset at the Attach File step

Acknowledgments

I thank American Journal Experts (AJE) for the final English language editing. I am grateful to the reviewers for their comments, which provided me with the opportunity to improve the article.
==== Refs
References

1 Woodruff R. A simple method for approximating the variance of a complicated estimate J. Am. Stat. Assoc. 26 1971 411 414
2 Hansen M. Hurwitz W. Madow W. Sample Survey Method and Theory. Volume II. Theory 1953 John Wiley & Sons New York, USA
3 Cochran W.G. Sampling Techniques third ed. 1977 John Wiley & Sons New York, USA
4 Särndal C.E. Swensson B. Wretman J.H. Model Assisted Survey Sampling 1992 Springer New York, USA
5 Tillé Y. Sampling and Estimation from Finite Populations 2020 John Wiley & Sons Hoboken, USA
6 Wolter K. Introduction to Variance Estimation second ed. 2007 Springer New York, USA
7 Aubry P. On the non-recursive implementation of multistage sampling without replacement MethodsX 8 2021 101553 34754820
8 Aubry P. Quaintenne G. Dupuy J. Francesiaz C. Guillemain M. Caizergues A. On using stratified two-stage sampling for large-scale multispecies surveys Ecol. Inform. 77 2023 102229
9 Sukhatme P.V. Sukhatme B.V. Sukhatme S. Asok C. Sampling Theory of Surveys with Applications third ed. 1984 Iowa State University Press Ames, Iowa, USA
10 Mukhopadhyay P. Theory and Methods of Survey Sampling second ed. 2009 PHI Learning New Delhi, India
11 Cornfield J. On samples from finite populations J. Am. Stat. Assoc. 39 1944 236 239
12 Aubry P. On univariate optimal partitioning by complete enumeration MethodsX 10 2023 102154 37091960
13 Raj D. Chandhok P. Sample Survey Theory 1998 Narosa Publishing House New Delhi, India
14 Konijn H. Statistical Theory of Sample Survey Design and Analysis 1973 North-Holland Amsterdam, The Netherlands
15 Lehtonen R. Pahkinen E. Practical Methods for Design and Analysis of Complex Surveys second ed. 2004 John Wiley & Sons Chichester, UK
16 Mukhopadhyay P. Complex Surveys. Analysis of Categorical Data 2016 Springer Singapore
17 Foreman E. Survey Sampling Principles 1991 Marcel Dekker New York, USA
18 Tepping B. Variance estimation in complex surveys Proceedings of the Social Statistics Section 1968 American Statistical Association (Alexandria, Virginia, USA) 11 18
19 Valliant R. Dever J. Kreuter F. Practical Tools for Designing and Weighting Survey Samples 2013 Springer New York, USA
20 SAS Institute.  SAS/STAT 15.3. User’s Guide 2023 SAS Institute Cary, North Carolina, USA
21 Carlson B. Software for sample survey data Armitage P. Colton T. Encyclopedia of Biostatistics second ed. 2005 John Wiley & Sons Chichester, UK
22 Hansen M. Hurwitz W. Madow W. Sample Survey Method and Theory. Volume I. Methods and Applications 1953 John Wiley & Sons New York, USA
23 Rust K. Kalton G. Strategies for collapsing strata for variance estimation J. Off. Stat. 3 1987 69 81
24 Zardetto D. ReGenesees: an advanced R system for calibration, estimation and sampling error assessment in complex sample surveys J. Off. Stat. 31 2015 177 203
25 Aubry P. Francesiaz C. Guillemain M. On the impact of preferential sampling on ecological status and trend assessment Ecol. Model. 492 2024 110707
26 Fisher W. On grouping for maximum homogeneity J. Am. Stat. Assoc. 53 1958 789 798
27 Hartigan J. Clustering Algorithms 1975 John Wiley & Sons New York, USA
28 Aubry P. On the correct implementation of the Hanurav–Vijayan selection procedure for unequal probability sampling without replacement Commun. Stat. Simul. Comput. 52 2023 1849 1877
29 Nijenhuis A. Wilf H. Combinatorial Algorithms for Computers and Calculators second ed. 1978 Academic Press New York, USA
30 Benedetti R. Espa G. Lafratta G. A tree-based approach to forming strata in multipurpose business surveys Survey Methodol. 34 2008 195 203
31 Hemelrijk J. Underlining random variables Stat. Neerl. 20 1966 1 7
32 Bohrnstedt G. Goldberger A. On the exact covariance of products of random variables J. Am. Stat. Assoc. 64 1969 1439 1442
33 Goodman L. On the exact variance of products J. Am. Stat. Assoc. 55 1960 708 713
