In many situations, the elements of a cluster may be too similar, so analyzing all the elements that make up the cluster would be a waste of resources. In these cases, it may be cheaper to select more clusters and take a subsample within each of them.
Lohr (2000)
In the previous chapter, the natural grouping of elements in the population was used to save financial and logistical costs when planning a cluster sampling strategy. However, the savings in operational terms are reflected in a high price to pay with respect to the statistical efficiency of the strategy. One possible solution to reduce the variance is to increase the cluster sample size, a solution that would increase operating costs.
To maintain a balance between financial costs and the advantages of the sampling strategy, it is possible to take advantage of homogeneity within clusters and, in this way, avoid carrying out a census within each selected cluster and instead select a subsample within the selected cluster. Because the structural behavior of the characteristic of interest inside the clusters is homogeneous, an estimate of the cluster total would have small variance. Of course, because there is no access to a sampling frame of elements, an enumeration must be carried out to build a sampling frame of elements in each and only in the selected clusters. Once the sampling frame of elements within the clusters is available, the subsamples of elements are selected. Bautista (1998) states that the basic principle of multistage sampling can be defined as the hierarchical process that carries out the following steps \(l\) times: - Construction of \(l\) sampling frames of units (clusters in the first \(l-1\) stages of the sampling design and elements in the last stage). - Application of a sampling design and selection of the samples (or subsamples) from each sampling frame.
Note that the concept of a sampling unit has been introduced, referring either to clusters of elements or to the elements themselves. If the sampling design has three stages, for example, if the goal is to obtain estimates about the behavior of students in a given city and no sampling frame of students is available, it is possible in a first stage to build a sampling frame of each and every school in the city and select a sample of schools using a certain sampling design. Once the schools are selected, in a second stage, a sampling frame of academic levels within schools (courses or classes) is built and a sample of levels is selected. Thus, in the third and final stage, a sampling frame of elements is built; that is, of students belonging to each selected level, and a sample of elements that will be observed and measured is drawn.
It is interesting to observe how the population, in its natural state, is subdivided through a “hierarchical” behavior, which in this particular case takes the following form:
\[ \underbrace{\textbf{City}}_{\text{Population $U$}} \Longrightarrow \underbrace{\textbf{Schools}}_{\text{PSU}} \Longrightarrow \underbrace{\textbf{Levels}}_{\text{SSU}} \Longrightarrow \underbrace{\textbf{Students}}_{\text{TSU}} \]
As notation, the first cluster subdivision of the original population is called the Primary Sampling Unit or PSU; the Secondary Sampling Unit or SSU is the sub-subdivision of the population, that is, the subdivision of the PSUs. The Tertiary Sampling Unit or TSU corresponds to the elements of the target population, which in this particular case are the students in the city.
The final sampling units are not always elements. Thus, it is possible to plan a two-stage cluster design, meaning that the secondary sampling units are clusters, or it is also possible to apply a four-stage design of elements, where the final sampling units are elements; for example, Bautista (1998) presents the following case:
\[ \underbrace{\textbf{City}}_{\text{Population $U$}} \Longrightarrow \underbrace{\textbf{Section}}_{\text{PSU}} \Longrightarrow \underbrace{\textbf{Block}}_{\text{SSU}} \Longrightarrow \underbrace{\textbf{Dwelling}}_{\text{TSU}} \Longrightarrow \underbrace{\textbf{Person}}_{\text{FSU}} \]
The basic principle of a multistage sampling strategy is to build estimates from the bottom up. But for the results of estimation based on the sampling design to be applicable, the following two assumptions must be satisfied: - Invariance: this suggests that the selection probability of a sample of sampling units (clusters or elements) does not depend on the sampling design of the previous stage. - Independence: interpreted as meaning that the subsampling of any sampling unit is carried out independently of the other sampling units, in the same stage or in higher or lower stages.
For the rest of the chapter, it is implicitly assumed that these properties are satisfied at each sampling stage of the strategy. If the assumptions are not satisfied, the reader may consult the section on multiphase sampling in the Advanced Topics chapter. To further establish the internal philosophy of multistage sampling, it is necessary to study the simplest of all sampling designs of this class: two-stage sampling.
Also called two-stage sampling by Mahalanobis (1946), this sampling design estimates the total of each cluster \(t_i\) through a subsample within the selected clusters of the population. In estimating the parameters of interest, two sources of variability are found, one in each stage. That is, there is variability due to the selection of primary sampling units or clusters and, of course, there is also variability due to the selection of a sample of elements, or secondary sampling units, within the selected clusters.
Suppose that the population of elements \(U\) is divided into \(N_I\) primary sampling units, which define a partition of the population, also called \textbf{clusters} and denoted as \(U_I=\{U_1,\ldots,U_{N_I}\}\). The \(i\)-th cluster \(U_i\), \(i=1,\dots,N_I\), has size \(N_i\). Särndal et al. (1992) provide a general framework for two-stage sampling, such that - A sample \(s_I\) of primary sampling units is selected from \(U_I\) according to a sampling design \(p_I(s_I)\). Note that \(S_I\) represents the random sample of clusters such that \(Pr(S_I=s_I)=p_I(s_I)\). - For each cluster \(U_i\), \(i=1,\dots,N_I\), selected in sample \(s_I\), a sample \(s_i\) of elements is selected according to a sampling design \(p_i(s_i)\). Note that \(S_i\) represents the random sample of elements such that \(Pr(S_i=s_i)=p_i(s_i)\).
This two-stage sampling design must satisfy the two properties of invariance and independence. Invariance means that the sampling designs \(p_i(s_i)\) of the second stage do not depend on the result in the first stage; that is, the sampling design must always be the same within each of the primary sampling units.
\[ Pr(S_i = s_i\left|\right. S_I=s_I) = Pr(S_i = s_i). \]
Note that the above implies that \(p_i(\cdot|s_I)=p_I(\cdot)\).
Independence means that the sample selection process in the second stage within each primary sampling unit does not depend on the selection processes used in the remaining primary sampling units. That is, subsampling in a particular primary sampling unit is independent of subsampling in other primary sampling units1; therefore, for each random sample \(S_I\) in the first stage, we have
\[ Pr\left(\bigcup_{i \in s_I}s_i|s_I\right)=\prod_{i \in s_I}Pr(s_i|s_I) \]
If the sampling design in the first stage is with replacement, then a cluster may appear more than once, and subsampling must be carried out as many times as that primary unit appears in the realized sample \(s_I\); this guarantees that the properties of independence and invariance are satisfied. In terms of support, it is also possible to speak of three types of support, namely:
\[ \begin{aligned} Q&=\bigcup_{r=1}^{\#Q_I} \bigcup_{i\in s_I^{(r)}} s_i, \ \ \ \ \text{with} \ \ s_i \in Q^i \notag \\ &=\left\{ \bigcup_{i\in s_I^{(r)}} s_i, \ \ \ \ \text{with} \ \ s_i \in Q^i, r=1,\ldots,\#Q_I\right\} \end{aligned} \]
where \(s_I^{(r)}\) denotes the r-th possible sample in the first stage and the cardinality of \(Q\) is given by
\[ \#Q=\prod_{i\in U_I}\#Q^i \]
And the sample of elements, or secondary sampling units, is given by
\[ S=\bigcup_{i\in S_I}S_i,\ \text{with} \ \ S_i\in Q^i \]
with random sample size given by
\[ n(S)=\sum_{i\in S_I}n_i \]
The definition of the supports at each stage and in general allows us to state that the two-stage sampling design is a genuine sampling design.
The two-stage sampling design satisfies - \(p(s)\geq0\) for every \(s\in Q\) - \(\sum_{s\in Q}p(s)=1\)
Proof.
First, we have
\[ \begin{aligned} p(s)&={Pr( \text{Select } s_I \text{ in stage one and select } \bigcup_{i\in s_I} s_i\ \text{ in stage two})}\\ &=p_I(s_I)\underbrace{Pr\left(\bigcup_{i\in s_I}s_i|s_I\right)}_{Independence}\\ &=p_I(s_I)\prod_{i\in s_I}\underbrace{Pr(s_i|s_I)}_{Invariance}\\ &=p_I(s_I)\prod_{i\in s_I}p_i(s_i) \end{aligned} \]
and it is clear that \(p(s)\geq0\). Now, to prove the second property, we have
\[ \begin{aligned} \sum_{s\in Q}p(s)&=\sum_{r=1}^{\#Q_I}\sum_{s_I^{(r)}}p(s)\\ &=\sum_{r=1}^{\#Q_I}\sum_{s_I^{(r)}}p_I(s_I^{(r)})\prod_{i\in s_I^{(r)}}p_i(s_i)\\ &=\sum_{r=1}^{\#Q_I}p_I(s_I^{(r)})\underbrace{\sum_{s_I^{(r)}}\prod_{i\in s_I^{(r)}}p_i(s_i)}_{=1}\\ &=\sum_{r=1}^{\#Q_I}p_I(s_I^{(r)})=1 \end{aligned} \]
where the equivalence to one of the second summation in the third equality is obtained by analogy with the proof of Result 5.1.1, where the stratified design was defined as a product.
To illustrate the preceding result, together with the integration of the concepts of supports in each of the stages, the following example was designed using a two-stage sampling design without replacement.
Our example population \(U_I\), given by
\(U_I\)={U_1, U_2, U_3}
Suppose that a sample \(s_I\) of primary sampling units of size \(n_I=2\) is selected using a sampling design without replacement such that
\[ p_I(s_I)=\begin{cases} 0.5, & \text{if } s_I=\{U_1, U_2\},\\ 0.4, & \text{if } s_I=\{U_1, U_3\},\\ 0.1, & \text{if } s_I=\{U_2, U_3\} \end{cases} \]
Now, suppose that within each selected primary unit, a single element is selected according to the following sampling designs
\[ \begin{aligned} p_1(S_1\left|\right. S_I)&=\begin{cases} 0.5, & \text{if } s_1=\{\text{Yves}\},\\ 0.5, & \text{if } s_1=\{\text{Ken}\} \end{cases}\\\\ p_2(S_2\left|\right. S_I)&=\begin{cases} 0.9, & \text{if } s_2=\{\text{Erik}\},\\ 0.1, & \text{if } s_2=\{\text{Sharon}\} \end{cases}\\\\ p_3(S_3\left|\right. S_I)&=\begin{cases} 1.0, & \text{if } s_3=\{\text{Leslie}\} \end{cases} \end{aligned} \]
That is, the final sample size is \(n=2\). And the support of the first stage is given by
\[ Q_I=\left\{\{U_1,U_2\},\{U_1,U_3\},\{U_2,U_3\}\right\}, \]
and the supports of the second stage are given by \(Q^1=\{\{\text{Yves}\},\{\text{Ken}\}\}\), \(Q^2=\{\{\text{Erick}\},\{\text{Sharon}\}\}\), and \(Q^3=\{\{\text{Leslie}\}\}\). Given the above, the support \(Q\) is given by
\[ Q=\left\{\bigcup_{i\in s_I^{(1)}}s_i,\bigcup_{i\in s_I^{(2)}}s_i,\bigcup_{i\in s_I^{(3)}}s_i\right\}, \]
where
\[ \bigcup_{i\in s_I^{(1)}}s_i=\left\{\{\text{Yves}, \text{Erick}\},\{\text{Yves},\text{Sharon}\},\{\text{Ken},\text{Erick}\},\{\text{Ken},\text{Sharon}\}\right\}, \]
\[ \bigcup_{i\in s_I^{(2)}}s_i=\left\{\{\text{Erick}, \text{Leslie}\},\{\text{Sharon},\text{Leslie}\}\right\}, \]
y
\[ \bigcup_{i\in s_I^{(3)}}s_i=\left\{\{\text{Yves}, \text{Leslie}\},\{\text{Ken},\text{Leslie}\}\right\}. \]
The probabilities \(\prod_{i\in s_I}p_i(s_i)\) and \(p_I(s_I)\) for all possible samples are as follows:
p(s_1) X p(s_2) p(s_I) p(s)
Yves Erick 0.5 X 0.9 0.5 0.225
Yves Sharon 0.5 X 0.1 0.5 0.025
Ken Erick 0.5 X 0.9 0.5 0.225
Ken Sharon 0.5 X 0.1 0.5 0.025
Erick Leslie 0.9 X 1.0 0.1 0.090
Sharon Leslie 0.1 X 1.0 0.1 0.010
Yves Leslie 0.5 X 1.0 0.4 0.200
Ken Leslie 0.5 X 1.0 0.4 0.200
Total 1.000It is observed that \(p(s)\) is a genuine sampling design. Note that within each possible first-stage sample, the sum of probabilities is equal to one. For example, for \(S_I=\{U_1,U_2\}\), the possible samples in the second stage correspond to \(\{\text{Yves, Erick}\}\), \(\{\text{Yves, Sharon}\}\), \(\{\text{Ken, Erick}\}\), and \(\{\text{Ken, Sharon}\}\) with probabilities 0.45, 0.05, 0.45, and 0.05, respectively, and the sum of these probabilities is equal to one.
The population parameters of interest can be written as: - The population total,
\[ t_y=\sum_{k \in U}y_k=\sum_{i=1}^{N_I}\sum_{k\in U_i}y_k=\sum_{i=1}^{N_I}t_{yi} \]
where \(t_{yi}=\sum_{k\in U_i}y_k\) is the total of the \(i\)-th primary sampling unit \(i=1,\dots,N_I\). - The population mean,
\[ \bar{y}_U=\frac{\sum_{k \in U}y_k}{N}=\frac{1}{N}\sum_{i=1}^{N_I}\sum_{k\in U_i}y_k=\frac{1}{N}\sum_{i=1}^{N_I}N_i\bar{y}_i \]
where \(\bar{y}_i=\dfrac{1}{N_i}\sum_{k\in U_i}y_k\) is the mean of the \(i\)-th primary sampling unit \(i=1,\dots,N_I\).
Our example population \(U_I\), given by
\(U_I=\{U_1, U_2, U_3\}\)
Suppose that a sample \(s_I\) of primary sampling units of size \(n_I=2\) is selected. The subsampling in the second stage is such that in each primary sampling unit selected in the first stage, only one element is selected, so that the sample size of elements is two. Define the support \(Q\) of elements if the sample selection is with replacement.
In the first stage, the first- and second-order inclusion probabilities of the primary sampling units induced by the sampling design \(p_I(s_I)\) are given by \(\pi_{Ii}\) and \(\pi_{Iij}\), respectively, with \(i,j\in U_I\). Therefore, we have
\[ \Delta_{Iij}=\begin{cases} \pi_{Iij}-\pi_{Ii}\pi_{Ij}, & \text{if } i,j \in U_I,\\ \pi_{Ii}(1-\pi_{Ii}), & \text{if } i=j \in U_I. \end{cases} \]
In the second stage, the first- and second-order inclusion probabilities of the elements in the \(i\)-th primary sampling unit, \(i\in S_I\), induced by the sampling design \(p_i(s_i)\) and conditional on \(U_i\) being selected in the first-stage sample, are given by \(\pi_{k|i}\) and \(\pi_{kl|i}\), respectively, for \(k,l\in U_i\), with \(\pi_{k|i}=Pr(k\in S_i|U_i\in S_I)\) and \(\pi_{kl|i}=Pr(k\in S_i,l\in S_i|U_i\in S_I)\). Therefore, we have
\[ \Delta_{kl|i}=\begin{cases} \pi_{kl|i}-\pi_{k|i}\pi_{l|i}, & \text{if } k\neq l,\\ \pi_{k|i}(1-\pi_{k|i}), & \text{if } k=l. \end{cases} \]
In general, from the definition of inclusion probability, we have the following result.
The first-order inclusion probability of the \(k\)-th element of \(U\) is given by
\[ \begin{aligned} \pi_{k}=Pr(k\in S)&=Pr(k\in S_i\ \text{and}\ i\in S_I)\notag\\ &=Pr(k\in S_i|i\in S_I)Pr(i\in S_I)=\pi_{k|i}\pi_{Ii} \end{aligned} \]
The second-order inclusion probability is given by
\[ \pi_{kl}=\begin{cases} \pi_{Ii}\pi_{k|i}, & \text{if } k=l\in U_i,\\ \pi_{Ii}\pi_{k|i}, & \text{if } k\neq l\in U_i,\\ \pi_{Iij}\pi_{k|i}\pi_{l|j}, & \text{if } k\in U,l\in U_j (i\neq l). \end{cases} \]
With the preceding result, we can use the general form of the Horvitz-Thompson estimator to find its particular expression and its variance under a two-stage sampling design (Särndal et al. 1992). However, to find a faster way to calculate the variance of the estimator, we need to use some well-known results from probability theory. These have been widely used in the field of sampling, but only later did Hansen et al. (1953) publish these results applied to sampling. In general, the goal is to express: - The expectation of a random variable as the expected value of conditional expectations. - The variance of a random variable as the sum of the variance of conditional expectations and the expectation of conditional variances.
Let \(U\) and \(H\) be random variables, then:
\[ E_1(U)=E_2(E_1(U|H)) \]
y, a su vez,
\[ Var_1(U)=E_2(Var_1(U|H))+Var_2(E_1(U|H)) \]
where subscript 1 denotes the expectation or variance induced by the distribution function of the random variable \(U\), and subscript 2 denotes the expectation or variance induced by the distribution function of the random variable \(H\).
Proof.
It is necessary to recall that \(Pr(U=U_i|H_j)=Pr(U=U_i,H=H_j)/Pr(H_j)\) and also that \(Pr(U=U_i)=\sum_j (U=U_i, H=H_j)\); consequently. - Expectation:
\[ \begin{aligned} E_1(U)&=\sum_i U_iPr(U=U_i)\\ &=\sum_i U_i\sum_j Pr(U=U_i, H=H_j)\\ &=\sum_i U_i\sum_j Pr(U=U_i|H=H_j)Pr(H=H_j)\\ &=\sum_j Pr(H=H_j) \sum_i U_iPr(U=U_i|H=H_j)\\ &=\sum_j Pr(H=H_j) E_2(U|H=H_j)\\ &=E_2(E_1(U|H)) \end{aligned} \] - Covariance: let W be a random variable and take \(x=E_2(U)\) and \(y=E_2(W)\)
\[ \begin{aligned} Cov(U,W)&=E(UW)-E(U)E(W)\\ &=E_1(E_2(UW))-E_1(E_2(U))E_1(E_2(W))\\ &=E_1(E_2(UW))-E_1(x)E_1(y)\\ &=E_1\left[E_2(UW)-xy\right]+E(xy)-E_1(x)E_1(y)\\ &=E_1\left[Cov_2(U,W)\right]+Cov_1(x,y)\\ &=E_1\left[Cov_2(U,W)\right]+Cov_1\left[E_2(U),E_2(W)\right] \end{aligned} \] - Variance: since variance is a particular case of covariance, then:
\[ \begin{aligned} Var(U)=Cov(U,U)&=E_1\left[Cov_2(U,U)\right]+Cov_1\left[E_2(U), E_2(U)\right]\\ &=E_1[Var_2(U)]+Var_1[E_2(U)] \end{aligned} \]
With the help of the preceding result, it is possible to obtain expressions for the Horvitz-Thompson estimator that show the variation in each of the two stages of this sampling design. The form taken by both the generic estimator and its corresponding variance is interesting because, since there are two sampling stages, the totals of the clusters are estimated in the first, and in the second stage the grand total is estimated using those estimates in the selected primary units. Since the estimation process is carried out in two stages, it is expected that there will be two sources of variation: the first due to estimating the totals of the primary sampling units and the second due to estimating the grand total. Assuming that four primary sampling units were selected, there would then be four estimates whose variance would be summarized in a single expression, while, on the other hand, there would be another source of variation when estimating the grand total.
Under two-stage sampling, the Horvitz-Thompson estimator is unbiased for the population total and takes the form
\[ \hat{t}_{y,\pi}=\sum_{i\in S_I}\sum_{k\in S_i}\frac{y_k}{\pi_{Ii}\pi_{k|i}}=\sum_{i\in S_I}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}} \]
with variance given by
\[ Var_{TS}(\hat{t}_{y,\pi})=\underbrace{\sum\sum_{U_I}\Delta_{Iij}\frac{t_i}{\pi_{Ii}}\frac{t_j}{\pi_{Ij}}}_{Var(PSU)}+\underbrace{\sum_{i\in U_I}\frac{Var_{p_i}(\hat{t_i})}{\pi_{Ii}}}_{Var(SSU)} \]
whose unbiased estimate is
\[ \widehat{Var}_{TS}(\hat{t}_{y,\pi})=\underbrace{\sum\sum_{S_I}\frac{\Delta_{Iij}}{\pi_{Iij}}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\frac{\hat{t}_{yj,\pi}}{\pi_{Ij}}}_{\widehat{Var}(PSU)} +\underbrace{\sum_{i\in S_I}\frac{\widehat{Var}(\hat{t}_{yi,\pi})}{\pi_{Ii}}}_{\widehat{Var}(SSU)} \]
where
\[ Var(\hat{t_i})=\sum\sum_{U_i}\Delta_{kl|i}\frac{y_k}{\pi_{k|i}}\frac{y_l}{\pi_{l|i}} \]
\[ \hat{t}_{yi,\pi}=\sum_{k\in S_i}\frac{y_k}{\pi_{k|i}} \]
representing the estimate of the total of the characteristic of interest in the \(i\)-th primary sampling unit, and
\[ \widehat{Var}(\hat{t_i})=\sum\sum_{S_i}\frac{\Delta_{kl|i}}{\pi_{kl|i}}\frac{y_k}{\pi_{k|i}}\frac{y_l}{\pi_{l|i}} \]
Note that the variation of the estimator is decomposed into the two stages specific to this design. It is also important to keep in mind that \(\widehat{Var}(PSU)\) and \(\widehat{Var}(SSU)\) are not unbiased estimators for \(Var(PSU)\) and \(Var(SSU)\), respectively. However, the entire expression \(\widehat{Var}_{TS}(\hat{t}_{y,\pi})\) is unbiased for \(Var_{TS}(\hat{t}_{y,\pi})\).
Proof.
To develop the preceding result, it is necessary to handle the two concepts inherent to sampling in two or more stages. a) Invariance: to select the primary sampling units, the same design must be used, and b) Independence: whatever design is chosen to select the elements within a primary sampling unit, it must not affect subsampling in any other primary sampling unit; therefore, any covariance existing at this stage will be null.
First, the Horvitz-Thompson estimator has the following form:
\[ \begin{aligned} \hat{t}_{y,\pi}&=\sum_{k\in S}\frac{y_k}{\pi_k}\\ &=\sum_{i\in S_I}\sum_{k\in S_i}\frac{y_k}{\pi_{Ii}\pi_{k|i}}\\ &=\sum_{i\in S_I}\frac{1}{\pi_{Ii}}\sum_{k\in S_i}\frac{y_k}{\pi_{k|i}}\\ &=\sum_{i\in S_I}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}} \end{aligned} \] - Unbiasedness of the estimator:
\[ \begin{aligned} E_p(\hat{t}_{y,\pi})&=E_{p_I}\left(E_p\left[\sum_{i\in S_I}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\left|\right. S_I\right]\right)\\ &=E_{p_I}\left(\sum_{i\in S_I}\underbrace{E_p\left[\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\left|\right. S_I\right]}_{\text{invariance}}\right)\\ &=E_{p_I}\left(\sum_{i\in S_I}\frac{E_{p_i}(\hat{t}_{yi,\pi})}{\pi_{Ii}}\right)\\ &=E_{p_I}\left(\sum_{i\in S_I}\frac{t_{yi,\pi}}{\pi_{Ii}}\right)\\ &=\sum_{i\in U_I}\frac{t_{yi,\pi}}{\pi_{Ii}}E_{P_I}(I_{Ii}(S_I))=t_y \end{aligned} \] - Variance:
\[ \begin{aligned} Var_p(\hat{t}_{y,\pi})=\underbrace{Var_{p_I}\left(E_p\left[\hat{t}_{y,\pi}\left|\right. S_I\right]\right)}_{Var(PSU)}+ \underbrace{E_{p_I}\left(Var_p\left[\hat{t}_{y,\pi}\left|\right. S_I\right]\right)}_{Var(SSU)} \end{aligned} \]
The first summand is equivalent to
\[ \begin{aligned} Var_{p_I}\left(E_p\left[\hat{t}_{y,\pi}\left|\right. S_I\right]\right)&=Var_{p_I}\left(E_p\left[\sum_{i\in S_I}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\left|\right. S_I\right]\right)\\ &=Var_{p_I}\left(\sum_{i\in S_I}\underbrace{\frac{E_p(\hat{t}_{y,\pi}\left|\right. S_I)}{\pi_{Ii}}}_{Invariance}\right)\\ &=Var_{p_I}\left(\sum_{i\in S_I}\frac{E_p(\hat{t}_{y,\pi})}{\pi_{Ii}}\right)\\ &=Var_{p_I}\left(\sum_{i\in S_I}\frac{t_{yi,\pi}}{\pi_{Ii}}\right)\\ &=\sum\sum_{U_I}\Delta_{Iij}\frac{t_{yi,\pi}}{\pi_{Ii}}\frac{t_{yj,\pi}}{\pi_{Ij}} \end{aligned} \]
The second summand takes the following form
\[ \begin{aligned} E_{p_I}\left(Var_p\left[\hat{t}_{y,\pi}\left|\right. S_I\right]\right) &=E_{p_I}\left(Var_p\left[\sum_{i\in S_I}\frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\left|\right. S_I\right]\right)\\ &=E_{p_I}\left(\sum_{i\in S_I}\frac{Var_p(\hat{t}_{yi,\pi}\left|\right. S_I)}{\pi_{Ii}^2}\right)\\ &=E\left(\sum_{i\in S_I}\left[\frac{Var(\hat{t}_{yi,\pi})}{\pi_{Ii}^2}\right]\right)\\ &=E_{p_I}\sum_{i\in U_I}\frac{I_{Ii}(S_I)}{\pi_{Ii}^2}Var_{p_i}(\hat{t}_{yi,\pi})\\ &=\sum_{i\in U_I}\left[\frac{Var(\hat{t}_{yi,\pi})}{\pi_{Ii}}\right] \end{aligned} \]
Then, the variance of the estimator is given by expression (7.1.15). - Estimated variance: to verify that \(\widehat{Var}_{TS}(\hat{t}_{y,\pi})\) is an unbiased estimator of the variance of the Horvitz-Thompson estimator, it must be kept in mind that
\[ \begin{aligned} E\left(\hat{t}_{yi,\pi}\hat{t}_{yj,\pi}\left|\right. S_I\right)&= \begin{cases} Var_{p_i}(\hat{y}_{yi,\pi})+(E_{p_i}(\hat{y}_{yi,\pi}))^2, & \text{if } i=j,\\ E_{p_i}(\hat{y}_{yi,\pi})E_{p_j}(\hat{y}_{yj,\pi}), & \text{if } i\neq j \end{cases} \notag\\ &= \begin{cases} Var(\hat{t}_{yi,\pi})+t_{yi,\pi}^2, & \text{if } i=j,\\ (t_{yi,\pi})(t_{yj,\pi}), & \text{if } i\neq j \end{cases} \end{aligned} \]
For the first part of the estimated variance, we have
\[ \begin{aligned} &\ \ E_{p_I}\left(E_p\left[\sum\sum_{S_I}\frac{\Delta_{Iij}}{\pi_{Iij}} \frac{\hat{t}_{yi,\pi}}{\pi_{Ii}}\frac{\hat{t}_{yj,\pi}}{\pi_{Ij}}\left|\right. S_I\right]\right)\\ &=E_{p_I}\sum\sum_{S_I}\frac{\Delta_{Iij}}{\pi_{Iij}}\frac{E_p(\hat{t}_{yi,\pi}\hat{t}_{yj,\pi}\left|\right. S_I)}{\pi_{Ii}\pi_{Ij}}\\ &=E\left(\sum_{i\in S_I}\sum_{j\neq i\in S_I}\frac{\Delta_{Iij}}{\pi_{Iij}}\frac{(t_{yi,\pi})}{\pi_{Ii}}\frac{(t_{yj,\pi})}{\pi_{Ij}} +\sum_{S_I}\frac{\Delta_{Iii}}{\pi_{Iii}}\frac{Var(\hat{t}_{yi,\pi})+t_{yi,\pi}^2}{\pi_{Ii}^2}\right)\\ &=E\left(\sum_{i\in S_I}\sum_{j\in S_I}\frac{\Delta_{Iij}}{\pi_{Iij}}\frac{(t_{yi,\pi})}{\pi_{Ii}}\frac{(t_{yj,\pi})}{\pi_{Ij}} +\sum_{S_I}\frac{Var(\hat{t}_{yi,\pi})}{\pi_{Ii}^2}(1-\pi_{Ii})\right)\\ &=\sum_{i\in U_I}\sum_{j\in U_I}\Delta_{Iij}\frac{(t_{yi,\pi})}{\pi_{Ii}}\frac{(t_{yj,\pi})}{\pi_{Ij}} -\sum_{U_I}Var(\hat{t}_{yi,\pi})\left(1-\frac{1}{\pi_{Ii}}\right) \end{aligned} \]
For the second part of the estimated variance, we have
\[ \begin{aligned} &E\left(E\left[\sum_{i\in S_I}\frac{\widehat{Var}(\hat{t}_{yi,\pi})}{\pi_{Ii}}\left|\right. S_I\right]\right)\\ &=E\left(\sum_{i\in S_I}\frac{Var(\hat{t}_{yi,\pi})}{\pi_{Ii}}\right)\\ &=\sum_{i\in U_I}Var(\hat{t}_{yi,\pi})\\ &=\sum_{U_I}\frac{Var(\hat{t}_{yi,\pi})}{\pi_{Ii}}+\sum_{U_I}Var(\hat{t}_{yi,\pi})\left(1-\frac{1}{\pi_{Ii}}\right) \end{aligned} \]
Adding these two quantities gives the result. Note that by themselves these quantities are not unbiased for their population counterparts; however, we have that:
\[ \begin{aligned} E\left[\widehat{Var}(PSU)\right]+E\left[\widehat{Var}(SSU)\right]=Var(\hat{t}_{y,\pi}) \end{aligned} \]
Regarding the form taken by the variance of the Horvitz-Thompson estimator, Särndal et al. (1992) state that: - It is convenient to estimate the two variance components \(Var(PSU)\) and \(Var(SSU)\) separately to get an idea of the contribution of variability in each stage. - If \(\pi_{k|i}=\pi_{kl|i}=1\) for every \(k,l \in U_i\) and for every \(U_i\in S_I\), then \(Var(SSU)=0\) and this design takes the form of a cluster design. - If \(\pi_{Ii}=\pi_{Iij}=1\) for every \(i,j=1,\ldots,N_I\), then this design becomes a stratified design.
Using the information from Example 7.1.1, verify, through a lexical-graphic exercise, the unbiasedness of the Horvitz-Thompson estimator.
In simple random cluster sampling, each and every element belonging to the clusters selected in sample \(s_I\) was measured. However, because in most situations clusters tend to be very similar in the structural behavior of the characteristic of interest, incorporating elements that do not bring new information would be considered a waste of economic and logistical resources. For this, it is more economical to take a larger sample of primary sampling units and carry out subsampling within each of them.
This sampling design assumes that the population is divided into \(N_{I}\) primary sampling units, from which a sample \(s_{I}\) of \(n_{I}\) units is selected using a simple random sampling design. The subsampling within each selected primary unit is also simple random. That is, for each selected primary sampling unit \(i\in s_{Ih}\) of size \(N_i\), a sample \(s_i\) of elements of size \(n_i\) is selected.
When selecting samples of primary and secondary units without replacement, the sampling algorithms given in Chapter 2 are used, so the following steps must be carried out: - Separate the population into \(N_I\) primary sampling units using the cluster sampling frame. - Select \(n_I\) clusters using any of the methods presented in Section 3.2.1; that is, the negative coordination method or the Fan-Muller-Rezucha method. - For each primary unit selected in the first-stage sample \(s_I\), select \(n_i\) elements, \(i\in S_I\), using any of the methods presented in Section 3.2.1.
When the sampling design is simple random in both stages, the following first- and second-order inclusion probabilities are obtained
\[ \begin{aligned} \pi_{Ii} &= \frac{n_I}{N_I} \\ \pi_{Iij} &= \frac{n_I(n_I-1)}{N_I(N_I-1)} \end{aligned} \]
respectively. On the other hand, the inclusion probability of an element or secondary sampling unit belonging to the \(i\)-th primary sampling unit \(i\in U_I\) is given by
\[ \begin{aligned} \pi_{k} &= \frac{n_I}{N_I}\frac{n_i}{N_i} \end{aligned} \]
Once the sample of primary units \(s_I\) is selected, a complete enumeration of the elements belonging to it is carried out to build a sampling frame that allows the selection of a subsample for the corresponding measurement of each and every element belonging to the selected subsample. In the cluster random sampling design, the estimator of the population total \(t_y\) was given by \(\hat{t}_{y,\pi}=\dfrac{N_i}{n_i}\sum_{i\in S_I}t_{yi}\) because the exact totals of each selected cluster were known through a census conducted within them. On the other hand, in SRS-SRS two-stage sampling, because not all elements of the selected primary units are measured, these totals \(t_{yi}\) must be estimated using the following expression
\[ \hat{t}_{yi,\pi}=\frac{N_i}{n_i}\sum_{k\in S_i}y_k=N_i\bar{y}_{U_i} \]
The following result leads to an estimate of the parameter of interest.
Under SRS-SRS two-stage sampling, the Horvitz-Thompson estimator is unbiased for the population total and takes the form
\[ \hat{t}_{y,\pi}=\frac{N_{I}}{n_{I}}\sum_{i\in S_{I}}\frac{N_i}{n_i}\sum_{k\in S_i}y_k \]
with variance given by
\[ \begin{aligned} Var_{MS}(\hat{t}_{y,\pi})=\frac{N_{I}^2}{n_{I}}\left(1-\frac{n_{I}}{N_{I}}\right)S^2_{t_{y}U_I}+ \frac{N_{I}}{n_{I}}\sum_{i\in U_{I}}\frac{N_i^2}{n_i}\left(1-\frac{n_i}{N_i}\right)S^2_{y_{U_i}} \end{aligned} \]
whose unbiased estimate is
\[ \begin{aligned} \widehat{Var}_{MS}(\hat{t}_{y,\pi})=\frac{N_{I}^2}{n_{I}}\left(1-\frac{n_{I}}{N_{I}}\right)S^2_{\hat{t}_{y}S_I}+ \frac{N_{I}}{n_{I}}\sum_{i\in S_{I}}\frac{N_i^2}{n_i}\left(1-\frac{n_i}{N_i}\right)S^2_{y_{S_i}} \end{aligned} \]
where \(S^2_{t_{y}U_I}\) is the population variance of the totals \(t_{yi}\), \(i\in U_I\), of each and every primary sampling unit, and \(S^2_{y_{U_i}}\) is the population variance among the elements within each primary sampling unit. Similarly, \(S^2_{\hat{t}_{y}s_I}\) and \(S^2_{y_{s_i}}\).
The first term in (7.2.6) refers to the variability due to the first stage of the sampling design, while the second summand refers to the additional variance due to subsampling in the primary sampling units. Lohr (2000) states that, as in the case of the cluster sampling design, if the primary sampling units have different sizes, then the variability of the estimator can be very large. If the sizes \(N_i\) of the clusters \(i\in U_I\) are very different from each other, the variance component will be large even if the structural behavior of the characteristic of interest is constant in each primary unit.
As we move forward in the programmatic development of this text, we find that although the principles of estimation are the same, survey design and the estimation of the parameters of interest become more complex. Lohr (2000) states that the best way to design a survey is to review it after it has concluded because, once the survey has ended, it is possible to evaluate the effect of the primary sampling units on the final estimate and, in this way, know where more logistical resources should be allocated to obtain better information. But even when knowledge of the population is acceptable, the question of sample size always arises. In particular, how many primary sampling units should be selected in the sample? And how many elements or secondary sampling units should be selected in the subsampling within the primary sampling units?
For example, particularly in area surveys, the larger the size of the primary sampling unit, the more variability can be expected within it. However, if the primary unit is very large, the benefits of financial and logistical savings may be lost.
The objective of a good sample is to collect the largest amount of information at the lowest economic and operational cost. Suppose that the population is divided into \(N_{I}\) primary sampling units, from which a sample \(s_{I}\) of \(n_{I}\) units is selected. Each primary sampling unit contains exactly \(N_i=M\) elements or secondary sampling units. The subsampling is such that a sample of exactly \(n_i=m\) secondary sampling units is selected. Therefore, the population and sample sizes are given by
\[ N=N_IM \ \ \ \ \ \text{and}\ \ \ \ \ n=n_Im \]
respectively. Thus, the estimator of \(t_y\) can be written as
\[ \begin{aligned} \hat{t}_{y,\pi}=\frac{N_{I}}{n_{I}}\frac{M}{m}\sum_{i\in S_{I}}\sum_{k\in S_i}y_k \end{aligned} \]
and its variance as
\[ \begin{aligned} Var_{MS}(\hat{t}_{y,\pi})= \frac{N_{I}^2}{n_{I}}\left(1-\frac{n_{I}}{N_{I}}\right)S^2_{t_{y}U_I}+ \frac{N_{I}^2M^2}{n_{I}m}\left(1-\frac{m}{M}\right)\bar{S}^2_{y_{U_i}} \end{aligned} \]
where \(\bar{S}^2_{y_{U_i}}=(1/N_I)\sum_{i\in U_I}S^2_{y_{U_i}}\).
Using the results of the decomposition of the sums of squares, the variance of the two-stage strategy (2SRS) takes the following form
\[ Var_{2SRS}(\hat{t}_{y,\pi})= \frac{N_I^2M}{n_I}\left[\frac{1}{N_I-1}(SST-SSW)+\left(\frac{M}{m}-1\right)\frac{SSW}{N_I(M-1)}\right] \]
whereas the variance of the simple random strategy, with population size equal to \(N=M\times N_I\) elements and sample size equal to \(n=m\times n_I\) elements, can be written as
\[ Var_{SRS}(\hat{t}_{y,\pi})=\frac{N_I^2}{n_I}\left(1-\frac{n_I}{N_I}\right)M\frac{SST}{MN_I-1} \]
To find the optimal values of \(n_I\) and \(m\) that will be used in the first and second sampling stages so that, given a cost function, the variance of the estimator is minimized2, we have the following result.
Considering the following cost function
\[ C=c_1n_I+c_2n_Im \]
where \(c_1\) is the cost of building the sampling frame in each primary unit selected in sample \(s_I\), and \(c_2\) is the cost of collecting information on the characteristic of interest for the elements or secondary units selected by subsampling. The optimal values of \(n_I\) and \(m\) that minimize the variance of the estimator given by expression (7.2.6), subject to the total survey cost given by (7.2.11), are
\[ n_I=\frac{C}{c_1+c_2m} \]
y
\[ m=M\bar{S}^2_{y_{U_i}}\sqrt{\frac{c_1/c_2}{S^2_{t_{y}U_I}-M\bar{S}^2_{y_{U_i}}}} \]
Proof.
The quantity to minimize is given in expression (7.2.10), which is subject to the constraint of the cost function (7.2.11). Using the method of Lagrange multipliers, we have
\[ \begin{aligned} \mathcal{L}(n_I,m,\lambda)=\frac{N_{I}^2}{n_{I}}\left(1-\frac{n_{I}}{N_{I}}\right)S^2_{t_{y}U_I}+ \frac{N_{I}^2M^2}{n_{I}m}\left(1-\frac{m}{M}\right)\bar{S}^2_{y_{U_i}}\\ +\lambda(c_1n_I+c_2n_Im-C) \end{aligned} \]
Setting the partial derivatives to zero gives
\[ \begin{aligned} \frac{\partial\mathcal{L}}{\partial n_I}&=-\frac{N_{I}^2M^2}{n_{I}^2}\left(\frac{1}{m}-\frac{1}{M}\right)\bar{S}^2_{y_{U_i}}- \frac{N_{I}^2}{n_{I}^2}S^2_{t_{y}U_I}+c_1\lambda+c_2m\lambda=0\\ \frac{\partial\mathcal{L}}{\partial m}&=-\frac{N_{I}^2M^2}{n_{I}^2m^2}\bar{S}^2_{y_{U_i}}+c_2n_I\lambda=0 \end{aligned} \]
From (7.2.15), we have
\[ \begin{aligned} n_I^2=-\dfrac{N_{I}^2M^2\left(\dfrac{1}{m}-\dfrac{1}{M}\right)\bar{S}^2_{y_{U_i}}+N_I^2S^2_{t_{y}U_I}}{c_1\lambda+c_2m\lambda} \end{aligned} \]
From (7.2.16), we have
\[ \begin{aligned} n_I^2=-\dfrac{N_{I}^2M^2\bar{S}^2_{y_{U_i}}}{c_2m^2\lambda} \end{aligned} \]
Equating the previous equations and solving for \(m\) gives the proof of the result.
If \(\bar{S}^2_{y_{U_i}}\), the variability of the characteristic of interest within the primary units, is large, then \(m\) will be large. It should be emphasized that the results are valid if the cost function is correct.
When the sampling strategy uses the Horvitz-Thompson estimator, we can use its general form to find its variance under any sampling design. The expression for the variance of the Horvitz-Thompson estimator under two-stage sampling is given by
\[ Var(\hat{t}_{\pi})=\sum\sum_{UI}\Delta_{Iij}\frac{t_j}{\pi_{Ij}}\frac{t_i}{\pi_{Ii}}+\sum_{UI}V_i/\pi_{Ii} \]
whose unbiased estimate is
\[ \widehat{Var}_1(\hat{t}_{\pi})=\sum\sum_{sI}\frac{\Delta_{Iij}}{\pi_{Ii}} \frac{\hat{t}_i}{\pi_{Ii}} \frac{\hat{t}_j}{\pi_{Ij}}+\sum_{sI}\hat{V}_i/\pi_{Ii} \]
The preceding expression involves calculating the variances of the variables within each cluster. In a large-scale survey, this can become very tedious, costly, and also very time-consuming. Särndal et al. (1992) (p. 139) give a possible solution to the problem: keeping the first part of the variance estimator as a general estimator of it. Thus, a simple but biased estimator is
\[ \widehat{Var}_2(\hat{t}_{\pi})= \sum\sum_{sI}\frac{\Delta_{Iij}}{\pi_{Ii}}\frac{\hat{t}_i}{\pi_{Ii}}\frac{\hat{t}_j}{\pi_{Ij}} \]
The preceding estimator overestimates the variance for the primary sampling units, but it also does so with (7.2.19). Another possible solution for estimating the variance of the Horvitz-Thompson estimator is to assume that sampling in the first stage was carried out with replacement. Thus, the biased variance estimate would be given by
\[ \widehat{Var}_3(\hat{t}_{\pi})=\frac{1}{n(n-1)}\sum_{i=1}^{n}\left(\frac{\hat{t}_i}{p_{Ii}}-\hat{t}_{\pi}\right)^2 \]
A special case of the preceding term is obtained by assuming that \(\pi_k=np_k\), and if sampling in the first stage was simple random, then \(p_k=\frac{1}{N}\). The variance estimator under the previous condition is
\[ \widehat{Var}(\hat{t}_{\pi})=\frac{N^2}{n(n-1)}\sum_{i=1}^{n}\left(\hat{t}_i-\frac{\sum_{i=1}^n\hat{t}_i}{n}\right)^2 =\frac{N^2}{n}S^2_{\hat{t}_i} \]
Srinath and Hidiroglou (1980) propose a fast method for estimating the variance of the Horvitz-Thompson estimator. It assumes that the selection method in the second stage is SRS and is invariant in the first stage (the sample in the first stage can be selected using any design); this implies that this variance estimator is unbiased and is given by
\[ \widehat{Var}_4(\hat{t}_{\pi})=-\frac{1}{2}\sum\sum_{sI}\frac{\Delta_{Iij}}{\pi_{Ii}}\left(\breve{t'}_i\breve{t'}_j\right) \]
where \(\breve{t'}_j=\frac{\hat{t}'_j}{\pi_{Ij}}\) and \(\hat{t}'_j=\frac{N_j}{n'_j}\sum_{s'_j} y_k\), where \(s'_j\) denotes a sample of \(n'_j\) elements. The rule for determining \(n'_j\) and obtaining the estimator \(\widehat{Var}_4\) is
\[ n'_j=\frac{n_i(1-\pi_{Ii})}{1-\pi_{Ii}(n_i/N_i)} \]
Simulation: the data from the 1996 FAMEX household expenditure survey (Canada Family Expenditure) were used. It has a total of 691 individuals and is divided into five clusters. The expenditure variable was used to estimate the total in a two-stage sample, and the FAMEX 1996 data, although they are survey data, were taken as the data of a universe.
The study aims to verify the results obtained above. For the sample design, three clusters were selected in the first stage; for each selected cluster, a sample whose size was 40% of that cluster was drawn. The sampling and subsampling were simple random SRS-SRS. The population total for the variable of interest is USD 711623, and the variance of the \(\pi\) estimator under the preceding conditions is 6595944566.
Thus, the following estimators were calculated for the variance of the estimated total \(\hat{t}_{\pi}\): - \(\widehat{Var}_1(\hat{t}_{\pi})\): the classical estimator when using two-stage sampling. - \(\widehat{Var}_2(\hat{t}_{\pi})\): corresponding to the first summand of the preceding estimator. - \(\widehat{Var}_3(\hat{t}_{\pi})\): the estimator assuming sampling with replacement. - \(\widehat{Var}_4(\hat{t}_{\pi})\): the estimator proposed by (Srinath and Hidiroglou 1980) (1.5).
The process was repeated \(B=5000\) times. The simulation was programmed in the R statistical package. In the simulation, the performance of an estimator \(\hat{V}\) was evaluated using its relative bias, \(RB\), and relative efficiency, \(RE\), defined as:
\[ RB=B^{-1}\sum_{b=1}^{B}\frac{\hat{V}_{b}-V}{V} \]
\[ RE=\frac{MSE(\hat{V}_{\pi})}{MSE(\hat{V})}, \]
where
\[ MSE(\hat{V})=B^{-1}\sum_{b=1}^{B}(\hat{V}_{b}-V)^2 \]
and \(\hat{V}_{b}\) was calculated in the b-th simulated sample. As can be seen, the classical estimator used under two-stage sampling, \(\hat{V}_{\pi}\), was used as the comparison baseline. Large values of RE(\(>1\)) represent high efficiency of the estimator \(\hat{V}\) compared with the classical estimator.
| \(\widehat{Var}_1(\hat{t}_{\pi})\) | \(\widehat{Var}_2(\hat{t}_{\pi})\) | \(\widehat{Var}_3(\hat{t}_{\pi})\) | \(\widehat{Var}_4(\hat{t}_{\pi})\) |
|---|---|---|---|
| 0.0008138860 | 0.2458789480 | -1.5021980054 | -0.0008792021 |
The empirical results indicate that the variance estimator for the Horvitz-Thompson estimator is unbiased, as is the estimator proposed by (Srinath and Hidiroglou 1980). However, estimators 2 and 3 have important relative bias, especially the one that assumes sampling with replacement; it can also be observed that the estimator of the first part of (7.2.20), although biased, has small magnitude. In particular, it is recommended to continue working with the classical estimator because computational advances now allow it. The relative efficiency of all estimators turned out to be negligible.
In the previous chapter, a cluster sampling design was carried out whose main characteristic is that the units within each cluster have relatively similar behavior. This led the estimates to be very far from reality because a sampling design that induced constant inclusion probabilities was used, even though the behavior of the cluster totals was not constant for the characteristics of interest.
On this occasion, we again face the difficulty of obtaining a sample of firms in the industrial sector without having a sampling frame that allows the direct inclusion of firms in the sample. However, it is possible to use as a base the area sampling proposed in the previous chapter, but the major difference is that, instead of a census in the selected geographic areas, we will carry out subsampling. Recall that the city is divided into five geographic zones labeled Zone A, located in the south; Zone B, located in the north; Zone C, located in the east; Zone D, located in the west; and Zone E, located in the center.
Suppose that there is no information about how many firms belong to each geographic zone, so it is not possible to carry out a self-weighting design. To guarantee good precision, it has been decided to select a simple random sample of four geographic zones, or primary sampling units. This is done using the sample function, although it is also acceptable to do it with the S.SI function from the TeachingSampling package.
data(BigLucy)
UI <- levels(BigLucy$Zone)
NI <- length(UI)
nI <- 20
samI <- S.SI(NI, nI)
selected_zones <- UI[samI]
selected_zones [1] "County13" "County18" "County19" "County20" "County25" "County37"
[7] "County38" "County44" "County50" "County51" "County52" "County60"
[13] "County68" "County7" "County76" "County80" "County83" "County84"
[19] "County9" "County92"Once the random draw is carried out, the selected geographic zones are: Zone B, Zone C, Zone D, and Zone E. The next step is the enumeration of each firm in the industrial sector belonging to each zone included in the sample. That is, a field operation must be planned in order to build a sampling frame for each primary unit. In total, four sampling frames of firms must be obtained.
Lucy1 <- BigLucy[which(BigLucy$Zone == selected_zones[1]), ]
Lucy2 <- BigLucy[which(BigLucy$Zone == selected_zones[2]), ]
Lucy3 <- BigLucy[which(BigLucy$Zone == selected_zones[3]), ]
Lucy4 <- BigLucy[which(BigLucy$Zone == selected_zones[4]), ]
Lucy5 <- BigLucy[which(BigLucy$Zone == selected_zones[5]), ]
Lucy6 <- BigLucy[which(BigLucy$Zone == selected_zones[6]), ]
Lucy7 <- BigLucy[which(BigLucy$Zone == selected_zones[7]), ]
Lucy8 <- BigLucy[which(BigLucy$Zone == selected_zones[8]), ]
Lucy9 <- BigLucy[which(BigLucy$Zone == selected_zones[9]), ]
Lucy10 <- BigLucy[which(BigLucy$Zone == selected_zones[10]), ]
Lucy11 <- BigLucy[which(BigLucy$Zone == selected_zones[11]), ]
Lucy12 <- BigLucy[which(BigLucy$Zone == selected_zones[12]), ]
Lucy13 <- BigLucy[which(BigLucy$Zone == selected_zones[13]), ]
Lucy14 <- BigLucy[which(BigLucy$Zone == selected_zones[14]), ]
Lucy15 <- BigLucy[which(BigLucy$Zone == selected_zones[15]), ]
Lucy16 <- BigLucy[which(BigLucy$Zone == selected_zones[16]), ]
Lucy17 <- BigLucy[which(BigLucy$Zone == selected_zones[17]), ]
Lucy18 <- BigLucy[which(BigLucy$Zone == selected_zones[18]), ]
Lucy19 <- BigLucy[which(BigLucy$Zone == selected_zones[19]), ]
Lucy20 <- BigLucy[which(BigLucy$Zone == selected_zones[20]), ]
LucyI <- rbind(
Lucy1, Lucy2, Lucy3, Lucy4, Lucy5, Lucy6, Lucy7, Lucy8, Lucy9,
Lucy10, Lucy11, Lucy12, Lucy13, Lucy14, Lucy15, Lucy16, Lucy17,
Lucy18, Lucy19, Lucy20
)
N1 <- dim(Lucy1)[1]
N2 <- dim(Lucy2)[1]
N3 <- dim(Lucy3)[1]
N4 <- dim(Lucy4)[1]
N5 <- dim(Lucy5)[1]
N6 <- dim(Lucy6)[1]
N7 <- dim(Lucy7)[1]
N8 <- dim(Lucy8)[1]
N9 <- dim(Lucy9)[1]
N10 <- dim(Lucy10)[1]
N11 <- dim(Lucy11)[1]
N12 <- dim(Lucy12)[1]
N13 <- dim(Lucy13)[1]
N14 <- dim(Lucy14)[1]
N15 <- dim(Lucy15)[1]
N16 <- dim(Lucy16)[1]
N17 <- dim(Lucy17)[1]
N18 <- dim(Lucy18)[1]
N19 <- dim(Lucy19)[1]
N20 <- dim(Lucy20)[1]
Ni <- c(
N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13, N14, N15,
N16, N17, N18, N19, N20
)
ni <- round(Ni * 0.12)
ni [1] 234 234 54 40 40 174 234 54 40 74 174 40 234 174 37 20 117 27 54
[20] 87sum(ni)[1] 2142When the first sampling stage concludes, it is known how many firms in the industrial sector belong to each geographic zone included in the sample. Zone B has 727 firms, Zone C has 974 firms, Zone D has 223 firms, and finally Zone E has a total of 165 firms. It has been decided that the sample sizes will correspond to a percentage of the size of each primary sampling unit. The sample size is 410 firms.
With the help of each of the four sampling frames, a simple random sample of firms is drawn according to the sizes established previously. Once the samples have been selected, they are combined using the rbind function, which simply merges the databases of the firms included in the sample.
sam1 <- sample(N1, ni[1])
sam2 <- sample(N2, ni[2])
sam3 <- sample(N3, ni[3])
sam4 <- sample(N4, ni[4])
sam5 <- sample(N5, ni[5])
sam6 <- sample(N6, ni[6])
sam7 <- sample(N7, ni[7])
sam8 <- sample(N8, ni[8])
sam9 <- sample(N9, ni[9])
sam10 <- sample(N10, ni[10])
sam11 <- sample(N11, ni[11])
sam12 <- sample(N12, ni[12])
sam13 <- sample(N13, ni[13])
sam14 <- sample(N14, ni[14])
sam15 <- sample(N15, ni[15])
sam16 <- sample(N16, ni[16])
sam17 <- sample(N17, ni[17])
sam18 <- sample(N18, ni[18])
sam19 <- sample(N19, ni[19])
sam20 <- sample(N20, ni[20])
sample_data1 <- Lucy1[sam1, ]
sample_data2 <- Lucy2[sam2, ]
sample_data3 <- Lucy3[sam3, ]
sample_data4 <- Lucy4[sam4, ]
sample_data5 <- Lucy5[sam5, ]
sample_data6 <- Lucy6[sam6, ]
sample_data7 <- Lucy7[sam7, ]
sample_data8 <- Lucy8[sam8, ]
sample_data9 <- Lucy9[sam9, ]
sample_data10 <- Lucy10[sam10, ]
sample_data11 <- Lucy11[sam11, ]
sample_data12 <- Lucy12[sam12, ]
sample_data13 <- Lucy13[sam13, ]
sample_data14 <- Lucy14[sam14, ]
sample_data15 <- Lucy15[sam15, ]
sample_data16 <- Lucy16[sam16, ]
sample_data17 <- Lucy17[sam17, ]
sample_data18 <- Lucy18[sam18, ]
sample_data19 <- Lucy19[sam19, ]
sample_data20 <- Lucy20[sam20, ]
sample_data <- rbind(
sample_data1, sample_data2, sample_data3, sample_data4, sample_data5, sample_data6,
sample_data7, sample_data8, sample_data9, sample_data10, sample_data11, sample_data12,
sample_data13, sample_data14, sample_data15, sample_data16, sample_data17, sample_data18,
sample_data19, sample_data20
)
head(sample_data) ID Ubication Level Zone Income Employees Taxes SPAM
52575 AB0000052575 C0144441K0157456 Medium County13 610 73 18.0 yes
53249 AB0000053249 C0283823K0018074 Small County13 77 81 0.5 yes
52071 AB0000052071 C0012349K0289548 Medium County13 990 145 47.0 no
53312 AB0000053312 C0159772K0142125 Small County13 175 44 1.0 yes
53109 AB0000053109 C0075639K0226258 Small County13 328 39 4.0 no
53511 AB0000053511 C0052610K0249287 Small County13 490 82 10.5 no
ISO Years Segments
52575 yes 38 County13 87
53249 no 34 County13 155
52071 yes 11 County13 37
53312 no 21 County13 161
53109 no 27 County13 141
53511 no 35 County13 181When information collection has concluded, the data file is loaded into the R environment and a data frame containing the values of the characteristics of interest in the general sample is constructed. In this particular case, it is called target_variables. Each firm included in the sample must carry the record indicating the geographic zone to which it belongs. For this exercise, the vector Area contains this information. Estimation in this two-stage sampling design is carried out using the E.2SI(NI,nI,Ni,ni,y,C) function, whose arguments are NI, the number of primary sampling units that make up the population; nI, the number of primary units included in sample \(s_I\); Ni, a vector of the sizes of the primary sampling units; ni, a vector containing the sample sizes in each primary sampling unit; y, the data file containing the information on the characteristics of interest; and finally C, a vector containing the membership of each secondary sampling unit in its corresponding primary unit.
target_variables <- data.frame(
Income = sample_data$Income,
Employees = sample_data$Employees,
Taxes = sample_data$Taxes
)
area <- as.factor(as.integer(sample_data$Zone))
E.2SI(NI, nI, Ni, ni, target_variables, area)The estimation results are shown in the following table. Note that with a similar sample size, the efficiency of this sampling strategy is much greater than that of a strategy using a cluster sampling design and is equivalent to that of a strategy using a simple random sampling design.
| N | Income | Employees | Taxes | |
|---|---|---|---|---|
| Estimation | 89110 | 35432820 | 5425504 | 893610 |
| Standard Error | 13532 | 5464183 | 822617 | 159249 |
| CVE | 15 | 15 | 15 | 18 |
| DEFF | Inf | 138 | 194 | 29 |
The gain in efficiency is due to the property of the two-stage design whereby, given an \(n\), it is possible to include more primary units in the first sampling stage. In this case, the number of clusters included in sample \(s_I\) is doubled, which decreases the variance component in the first stage. The variability component that dominates the variance in this estimate is the dispersion within the primary units and is due to the heterogeneity of the clusters.
The theory discussed so far in the previous sections is applicable when the primary sampling units are selected from a stratum. As will be seen later, there are no new estimation or design principles involved in the development of this sampling strategy when the goal is to estimate the total of the characteristic of interest \(t_y\) for a population divided into \(H\) strata.
It is assumed that sampling in each stratum respects the principle of independence. The estimates of the total, as well as the calculation and estimation of the variance, are simply the result of adding or summing the respective quantity for each stratum.
For example, suppose that within each stratum \(U_h\), \(h=1,\ldots, H\), there are \(N_{Ih}\) primary sampling units, from which a sample \(s_{Ih}\) of \(n_{Ih}\) units is selected using a simple random sampling design. Suppose also that the subsampling within each selected primary unit is also simple random. That is, for each selected primary sampling unit \(i\in s_{Ih}\) of size \(N_i\), a sample \(s_i\) of elements of size \(n_i\) is selected. When the secondary sampling units or elements are selected, the measurement process and the estimation process are carried out, for which the estimator of the total is given by the following result.
Under stratified SRS-SRS two-stage sampling, the Horvitz-Thompson estimator is unbiased for the population total and takes the form
\[ \hat{t}_{y,\pi}=\sum_{h=1}^H\hat{t}_{yh,\pi}=\sum_{h=1}^H\left[\frac{N_{Ih}}{n_{Ih}}\sum_{i\in S_{Ih}}\frac{N_i}{n_i}\sum_{k\in S_i}y_k\right] \]
with variance given by
\[ \begin{aligned} Var_{STMS}(\hat{t}_{y,\pi})&=\sum_{h=1}^HVar(\hat{t}_{yh,\pi})\\ &=\sum_{h=1}^H\left[\frac{N_{Ih}^2}{n_{Ih}}\left(1-\frac{n_{Ih}}{N_{Ih}}\right)S^2_{t_{yh}U_I}+ \frac{N_{Ih}}{n_{Ih}}\sum_{i\in S_{Ih}}\frac{N_i^2}{n_i}\left(1-\frac{n_i}{N_i}\right)S^2_{y_{U_i}}\right] \end{aligned} \]
whose unbiased estimate is
\[ \begin{aligned} \widehat{Var}_{STMS}(\hat{t}_{y,\pi})&=\sum_{h=1}^H\widehat{Var}(\hat{t}_{yh,\pi})\\ &=\sum_{h=1}^H\left[\frac{N_{Ih}^2}{n_{Ih}}\left(1-\frac{n_{Ih}}{N_{Ih}}\right)S^2_{\hat{t}_{yh}S_I}+ \frac{N_{Ih}}{n_{Ih}}\sum_{i\in S_{Ih}}\frac{N_i^2}{n_i}\left(1-\frac{n_i}{N_i}\right)S^2_{y_{S_i}}\right] \end{aligned} \]
where \(S^2_{t_{yh}U_I}\) is the population variance of the totals \(t_{yi}\), \(i\in U_I\), of each and every primary sampling unit within stratum \(h\), and \(S^2_{y_{U_i}}\) is the population variance among the elements within each primary sampling unit in stratum \(h\). Similarly, \(S^2_{\hat{t}_{yh}s_I}\) and \(S^2_{y_{s_i}}\).
This sampling design is used to improve the efficiency of the SRS-SRS strategy. Särndal et al. (1992) state that it is possible to stratify the population according to a size measure, so that sampling units with similar behavior are grouped in the same stratum. It is very interesting to note that a particular choice within the subsampling of the primary units would make the Horvitz-Thompson estimator very convenient to calculate. In fact, if for each primary unit \(i\in S_{I_h}\) selected in the sample of each stratum \(h\), \(h=1,\ldots,H\), we have
\[ c=\frac{n_i}{N_i}\frac{n_{Ih}}{N_{Ih}} \]
then the estimator takes the following form
\[ \hat{t}_{y,\pi}=\frac{1}{c}\sum_{h=1}^H\sum_{i\in S_{Ih}}\sum_{k\in S_i}y_{hik} \]
This means that, in the computational calculation of the estimate, the values of the characteristic of interest are simply summed regardless of the primary unit or stratum to which they belong. This class of estimators is known as self-weighting estimators. The quantity \(c\) has a very simple interpretation: it is the expected sampling fraction for the elements. Thus, if the goal is to select a sample with an average of 1% of secondary sampling units or elements selected in each stratum, then \(k=\dfrac{1}{100}\).
In many two-stage surveys, self-weighting designs are commonly found. This class of designs assumes that in the first sampling stage a sample \(S_I\) of primary sampling units is selected whose inclusion probabilities are proportional to their size, so that if \(N\) is the size of population \(U\) of secondary sampling units or elements and \(n\) is the resulting sample size, then
\[ \pi_{Ii}=\frac{N_i}{N}n_I \ \ \ \ \ i\in U_I \]
Later, in the second sampling stage, samples \(s_i\), \(i\in S_I\), of secondary units or elements of constant size \(n_i=n_0\) are selected for each primary unit included in the sample. Therefore, the inclusion probability of the secondary units will be
\[ \pi_{k|i}=\frac{n_0}{N_i} \ \ \ \ \ i\in S_I \]
Thus, the overall inclusion probability of the \(k\)-th element is constant and is given by
\[ \pi_k=\pi_{Ii}\pi_{k|i}=n_I\frac{N_i}{N}\frac{n_0}{N_i}=n_I\frac{n_0}{N}=\frac{n}{N}=c \ \ \ \ \ k\in U_i \]
and the Horvitz-Thompson estimator takes the following form
\[ \hat{t}_{y,\pi}=\sum_{k\in S}\frac{y_k}{\pi_k}=\frac{1}{c}\sum_{i\in S_I}\sum_{k\in S_i}y_k=\frac{N}{n}\sum_{k\in S}y_k \]
Note the ease of calculating the estimator. This class of self-weighting designs is used when fieldwork needs to be controlled, so the number of interviews in each primary unit included in the sample will be constant.
Särndal et al. (1992) state that despite their complexity, designs with three or more stages are widely used in large surveys. Two-stage sampling can be generalized through the following result, where it is assumed that there are \(r\) sampling stages. In this way, the population is divided into \(N_{I}\) primary sampling units, from which a sample \(s_{I}\) of \(n_{I}\) units is selected using a sampling design \(p_I(S_I)\). It is assumed that it is possible to construct an estimator3 \(\hat{t}_{yi}\) for each total \(t_{yi}\), \(i\in S_I\), of the selected primary units and that this estimator is unbiased for the remaining \(r-1\) stages of the sampling design. Therefore
\[ E(\hat{t}_{yi}\left|\right. S_I)=t_{yi} \]
Note that the final sampling units do not necessarily have to be elements; they may also be clusters. The principles of independence and invariance continue to hold in all stages of the sampling design. Thus, the foundation of this sampling design is the accumulation of estimates from the last stage up to the first. This is summarized in the following results in the next section.
Under \(r\)-stage sampling, the Horvitz-Thompson estimator is unbiased for the population total and takes the form
\[ \hat{t}_{y,\pi}=\sum_{i\in S_I}\frac{\hat{t}_{yi}}{\pi_{Ii}} \]
with variance given by
\[ Var_{TS}(\hat{t}_{y,\pi})=\underbrace{\sum\sum_{U_I}\Delta_{Iij}\frac{t_i}{\pi_{Ii}}\frac{t_j}{\pi_{Ij}}}_{Var(PSU)} +\underbrace{\sum_{i\in U_I}\frac{V_i}{\pi_{Ii}}}_{Var(Rest)} \]
whose unbiased estimate is
\[ \widehat{Var}_{TS}(\hat{t}_{y,\pi})=\underbrace{\sum\sum_{S_I}\frac{\Delta_{Iij}}{\pi_{Iij}}\frac{\hat{t}_{yi}}{\pi_{Ii}}\frac{\hat{t}_{yj}}{\pi_{Ij}}} _{\widehat{Var}(PSU)}+\underbrace{\sum_{i\in S_I}\frac{\hat{V}_i}{\pi_{Ii}}}_{\widehat{Var}(Rest)} \]
where \(V_i=Var(\hat{t}_{yi}\left|\right. S_I)\) and \(\hat{V}_i\) is an unbiased estimator of \(V_i\) such that \(E(\hat{V}_i\left|\right. S_I)=V_i\) for every \(i\in U_I\).
Proof.
This proof is carried out recursively by writing the estimator and the variance as a function of the unbiased estimators of the subsequent stages at lower levels. It must be kept in mind that Result 7.2.2 extends naturally. For example, for the three-stage design, we have
\[ Var(U)=V_1[E_2(E_3(U))]+E_1[V_2(E_3(U))]+E_1[E_2(V_3(U))] \]
A scheme used in practice because of the simplicity of the estimation process consists of selecting a sample of \(m_I\) primary sampling units using a sampling design with replacement that induces selection probabilities \(p_{Ii}\) with \(i\in U_I\) such that \(\sum_{i=1}^{N_i}p_{Ii}=1\). Within each primary sampling unit selected in the random draw with replacement, a subsample is taken (with or without replacement). Although there is a loss of efficiency when sampling is with replacement, it is offset by a logistical gain in the process of estimating the variances required for each characteristic of interest. According to Särndal et al. (1992), the general process of sampling with replacement is as follows: - In the first stage, a random sample is selected according to a sampling design with replacement such that \(p_{Ii}\) with \(i\in U_I\) is the selection probability of the \(i\)-th primary sampling unit. - In the following stages4, the properties of independence and invariance are maintained regardless of whether the design within the selected primary units is with or without replacement. - If a sampling unit is selected more than once, as many subsamples must be carried out as times it was selected in the first stage.
Under a multistage sampling design, the Hansen-Hurwitz estimator for the total \(t_{y}\), its variance, and its estimated variance are given by
\[ \hat{t}_{y,p}=\frac{1}{m_I}\sum_{v=1}^{m_I}\frac{\hat{t}_{{yi}_v}}{p_{Ii_v}} \]
\[ Var(\hat{t}_{y,p})=\frac{1}{m_I}\sum_{i=1}^{N_I}p_{Ii}\left(\frac{t_{yi}}{p_{Ii}}-t_y\right)^2+\frac{1}{m_I}\sum_{i=1}^{N_I}\frac{V_i}{p_{Ii}} \]
\[ \widehat{Var}(\hat{t}_{y,p})=\frac{1}{m_I(m_I-1)}\sum_{v=1}^{m_I}\left(\frac{\hat{t}_{yi_v}}{p_{Ii_v}}-\hat{t}_{y,p}\right)^2 \]
respectively. Here \(\hat{t}_{yi}\) is an unbiased estimator of the total of the characteristic of interest \(y\) in primary unit \(U_i\), \(i\in S_I\), and \(V_i=Var(\hat{t}_{yi}\left|\right. S_I)\) is the variance of \(\hat{t}_{yi}\) in the second stage. Note that \(\hat{t}_{y,p}\) is unbiased for \(t_y\) and that \(\widehat{Var}(\hat{t}_{y,p})\) is unbiased for \(Var(\hat{t}_{y,p})\).
Proof.
The proof begins by defining the random variables
\[ Z_v=t_{yi}/p_{Ii} \ \ \ \ \ \ i\in U_I \ \ \ v=1,\ldots,m_I \]
y
\[ \hat{Z}_v=\hat{t}_{yi}/p_{Ii} \ \ \ \ \ \ i\in U_I \ \ \ v=1,\ldots,m_I \]
Both \(Z_v\) and \(\hat{Z}_v\) are sequences of independent and identically distributed random variables. However, respecting the principles of independence and invariance, the expectation is given by
\[ \begin{aligned} E(\hat{Z}_v)=E(E(\hat{Z}_v\left|\right. S_I))=E(Z_v)=t_y \end{aligned} \]
and the variance is
\[ \begin{aligned} Var(\hat{Z}_v)&=Var(E(\hat{Z}_v\left|\right. S_I))+E(Var(\hat{Z}_v\left|\right. S_I))\\ &=Var(Z_v)+E(Var(\hat{t}_{yi}/p_{Ii}\left|\right. S_I))\\ &=Var(Z_v)+E(V_i/p_{Ii}^2)\\ &=\sum_{i=1}^{N_I}p_{Ii}\left(\frac{t_{yi}}{p_{Ii}}-t_y\right)^2+\sum_{i=1}^{N_I}\frac{V_i}{p_{Ii}} \end{aligned} \]
Now, since \(\hat{t}_{y,p}=\overline{\hat{Z}}\) and using Result 2.2.11, the unbiased estimator of the variance corresponds to the expression given in (7.4.8).
Given the simplification in calculating the variance, Bautista (1998) proposes using it even when the sampling design is without replacement. However, he warns that this estimator generally overestimates the variance, which leads to more conservative confidence intervals and slightly higher coefficients of variation.
| Municipality | \(N_i\) | \(n_i\) | \(y_k\) |
|---|---|---|---|
| 1 | 20 | 4 | 105 |
| 1 | 20 | 4 | 118 |
| 1 | 20 | 4 | 102 |
| 1 | 20 | 4 | 110 |
| 2 | 23 | 5 | 108 |
| 2 | 23 | 5 | 117 |
| 2 | 23 | 5 | 134 |
| 2 | 23 | 5 | 108 |
| 2 | 23 | 5 | 119 |
| 3 | 18 | 4 | 201 |
| 3 | 18 | 4 | 201 |
| 3 | 18 | 4 | 210 |
| 3 | 18 | 4 | 206 |
| 4 | 28 | 6 | 157 |
| 4 | 28 | 6 | 141 |
| 4 | 28 | 6 | 129 |
| 4 | 28 | 6 | 170 |
| 4 | 28 | 6 | 104 |
| 4 | 28 | 6 | 110 |
% ——————————————————————————–