Variance Component Estimation in Multistage Sampling

Special case: srswor at first and second stages

Suppose the first stage is an srswor of \(\small m\) out of \(\small M\) PSUs and the second stage is a sample of \(\small n_{i}\) elements selected by srswor from the population of \(\small N_{i}\). The \(\small \pi\)-estimator is

\[ \small \hat{t}_{\pi} = \frac{M}{m} \sum_{i \in s} \frac{N_{i} }{n_{i} } \sum_{k \in s_{i} } y_{k} \]

Its variance is equal to \[\begin{equation} \label{srs.relvar.gen} \small V\left(\hat{t}_{\pi } \right)=\frac{M^{2} }{m} \frac{M-m}{M} S_{U1}^{2} +\frac{M}{m} \sum _{i\in U}\frac{N_{i}^{2} }{n_{i} } \frac{N_{i} -n_{i} }{N_{i} } S_{U2i}^{2} \end{equation}\] where \(\small S_{U1}^{2} =\frac{\sum _{i\in U}\left(t_{i} -\bar{t}_{U} \right)^{2}}{M-1}\) with \(t_{i}\) being the population total of \(\small y\) in PSU \(\small i\) and \(\small \bar{t}_{U} ={\sum _{i\in U}t_{i} / M}\) is the mean total per PSU.

If \(\bar{n}\) elements are selected in each PSU and the sampling fractions of PSUs and elements within PSUs are all small, then the relvariance can be written as \[\begin{equation} \label{eq:relvarsrs} \small \frac{V(\hat{t}_{\pi})}{t_{U}^2} = \frac{B^2}{m} + \frac{W^2}{m\bar{n}} \end{equation}\]

where \(\small B^{2} ={S_{U1}^{2} / \bar{t}_{U}^{2} = M^{2}S_{U1}^2/t_{U}^2 }\) is the unit relvariance among PSU totals and \(\small W^{2} = M \sum _{i\in U}N_{i}^{2} S_{U2i}^{2}/t_{U}^2\). The term \(\small B^{2}\) is called the “between (PSU) component” while \(\small W^{2}\) is the “within component”. Expression \(\small \eqref{eq:relvarsrs}\) is the form used in the R function, BW2stageSRS. Textbooks often list a specialized form of \(\small \eqref{eq:relvarsrs}\) that requires that all PSUs have the same size, \(\small N_i \equiv \bar{N}\), and that \(\small \bar{n}\) elements are selected in each. In that case, the second-stage sampling fraction is \({\bar{n} / \bar{N}}\). This implies that the sample is self-weighting: \(\small \pi _{i} \pi _{k\left|i\right. } = m\bar{n}/M\bar{N}\). The relvariance based on \(\small \eqref{srs.relvar.gen}\) then simplifies to the less general form \[\begin{equation} \notag \small \frac{V(\hat{t}_{\pi})}{t_{U}^2} = \frac{1}{m} \frac{M-m}{M} B^{2} +\frac{1}{m\bar{n}} \frac{\bar{N}-\bar{n}}{\bar{N}} W^{2} \notag \end{equation}\]

where \(\small W^{2} =\frac{1}{M\bar{y}_{U}^{2} } \sum _{i\in U}S_{U2i}^{2}\).

Assuming that \(\small \bar{n}\) elements are selected in each sample PSU, and \(\small m/M\) and \(\small \bar{n}/N_{i}\) are both small, the more general form of the relvariance in \(\small \eqref{eq:relvarsrs}\) can also be written in terms of a measure of homogeneity \(\small \delta\) as follows:

\[\begin{equation} \label{eq:deltavar.2st} \small \frac{V\left(\hat{t}_{\pi } \right)}{t_{U}^{2} } \doteq \frac{\tilde{V}}{m\bar{n}} k \left[1+\delta \left(\bar{n}-1\right)\right] \end{equation}\] where \(\small \tilde{V}= S_{U}^2/\bar{y}_{U}^{2}\), \(\small k=(B^2 + W^2)/\tilde{V}\), and

\[\begin{equation} \label{eq:delta} \small \delta = \frac{B^{2} }{B^{2} +W^{2} }. \end{equation}\] With some effort, it can be shown that when \(\small N_{i}=\bar{N}\) and both \(\small M\) and \(\small \bar{N}\) are large,

\[\begin{equation} \small \frac{S_{U}^{2} }{\bar{y}_{U}^{2} } = {\frac{1}{\bar{y}_{U}^{2} } \frac{\sum _{i\in U}\sum _{k\in U_{i}}\left(y_{k} -\bar{y}_{U} \right)^{2} }{\left(N-1\right)} } \notag \\ \small {\doteq B^{2} +W^{2} } \notag \end{equation}\] i.e., the population relvariance can be written as the sum of between and within relvariances. If \(\small k=1\), \(\small \eqref{eq:deltavar.2st}\) equals the expression found in many textbooks. However, when the population count of elements per cluster varies, \(\small k\) may be far from 1, as will be illustrated in an example below. In those cases, \(\small \eqref{eq:deltavar.2st}\) with an estimate of the actual \(\small k\) should be used for determining sample sizes and computing advance estimates of coefficients of variation.

Expressions \(\small \eqref{eq:relvarsrs}\) and \(\small \eqref{eq:deltavar.2st}\) are useful for sample size calculation since the number of sample PSUs and sample units per PSU are explicit in the formula. Equation \(\small \eqref{eq:deltavar.2st}\) also connects the variance of the estimated total to the variance that would be obtained from a simple random sample since \(\small \tilde{V}/m\bar{n}\) is the relvariance of the estimated total in an srswor of size \(\small m\bar{n}\) when the sampling fraction is small. The product \(\small k[1 + \delta(\bar{n}-1)]\) is a type of design effect. When \(\small k=1\), the term \(\small 1+\delta \left(\bar{n}-1\right)\) is the approximate design effect found in many textbooks.

The next example uses the MDarea.popA from PracTools. This dataset is based on the U.S. Census counts from the year 2000 for Anne Arundel County in the US state of Maryland. The geographic divisions used in this dataset are called tracts and block groups. Tracts are constructed by the US Census Bureau to have a desired population size of 4,000 people. Block groups (BGs) are smaller with a target size of 1,500 people. Counts of persons in the dataset are the same for most tracts and block groups as in the 2000 Census.

• Example. Between and within variance components in srs/srs design The R function BW2stageSRS will calculate the unit relvariance of a population, \(\small B^{2} + W^{2}\) for comparison, the ratio \(\small k=(B^2+W^2)/(S_{U}^{2}/\bar{y}_{U}^2)\), and the full version of \(\small \delta\) in \(\small \eqref{eq:delta}\). The function assumes that the entire sampling frame is an input. The full R code for this example is in the file Example 9.2.R, available at Valliant, Dever, and Kreuter (2018a). We first compute the results using the PSU and SSU variables as clusters. These fields are created so that all PSUs have the same size; likewise, all SSUs have the same size. For the variable y1 in the Maryland population, the code is

 require(PracTools)
 data(MDarea.popA)
 BW2stageSRS(MDarea.popA$y1, psuID=MDarea.popA$PSU)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#> 0.007979563 1.455575239 1.463117997 1.463554802 1.000298544 0.005452179
 BW2stageSRS(MDarea.popA$y1, psuID=MDarea.popA$SSU)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#>  0.03661379  1.42824599  1.46311800  1.46485978  1.00119046  0.02499474

The values of \(\small \delta\) are 0.005 for PSU and 0.025 for SSU. Next, to illustrate the dramatic effect that varying sizes of clusters can have, we compute the same statistics as above using tracts and block groups (BGs) within tracts as clusters. These vary substantially in the number of persons in each cluster. A new variable called trtBG is computed since the values of the variable, BLKGROUP, are nested within each tract:

   trtBG <- 10*MDarea.popA$TRACT + MDarea.popA$BLKGROUP
   BW2stageSRS(MDarea.popA$y1, psuID=MDarea.popA$TRACT)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#>   0.2603356   1.8405833   1.4631180   2.1009189   1.4359190   0.1239151
   BW2stageSRS(MDarea.popA$y1, psuID=trtBG)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#>   0.3495410   1.9488250   1.4631180   2.2983660   1.5708685   0.1520824

The value of \(\small \delta\) is 0.124 TRACTs are clusters and 0.152 when trtBG defines clusters. The measures of homogeneity increase substantially when tracts or BGs are clusters compared to the PSU and SSU results. This is entirely due to the increase in \(\small B^{2}\) when units with highly variable sizes are used and an srs is selected. For example, \(\small B^{2} =0.0079\) for y1 when PSU is a cluster but is 0.2605 when TRACT is a cluster.

More General Two-stage Designs

Variances of estimators in two-stage designs more complicated than simple random sampling at each stage can be written as a sum of components. However, these have limited usefulness in determining sample sizes for the same reason that \(\small \eqref{eq:pivar}\) is not. A more convenient formulation is the case where PSUs are selected with varying probabilities but with replacement, and the sample within each PSU is selected by srswor. With-replacement designs may not often be used in practice but have simple variance formulae. The pwr-estimator of a total (Särndal, Swensson, and Wretman 1992) is \[ \small \hat{t}_{pwr} =\frac{1}{m} \sum _{i\in s}\frac{\hat{t}_{i} }{p_{i} } \] where \(\small \hat{t}_{i} =\frac{N_{i} }{n_{i} } \sum _{k\in s_{i} }y_{k}\) is the estimated total for PSU \(\small i\) from a simple random sample and \(\small p_{i}\) is the one-draw selection probability of PSU \(\small i\). The variance of \(\small \hat{t}_{pwr}\) is

\[\begin{equation} \label{eq:v2st} \small V\left(\hat{t}_{pwr} \right)=\frac{1}{m} \sum _{i\in U}p_{i} \left(\frac{t_{i} }{p_{i} } -t_{U} \right)^{2} +\sum _{i\in U}\frac{N_{i} ^{2} }{mp_{i} n_{i} } \left(1-\frac{n_{i} }{N_{i} } \right)S_{U2i}^{2}. \end{equation}\]

Making the assumption that \(\small \bar{n}\) elements are selected in each PSU, the variance reduces to \[ \small V\left(\hat{t}_{pwr} \right)=\frac{S_{U1\left(pwr\right)}^{2} }{m} +\frac{1}{m\bar{n}} \sum _{i\in U}\left(1-\frac{\bar{n}}{N_{i} } \right)\frac{N_{i}^{2} S_{U2i}^{2} }{p_{i} } \] where, in this case, \(S_{U1\left(pwr\right)}^{2} =\sum _{i\in U}p_{i} \left(\frac{t_{i} }{p_{i} } -t_{U} \right)^{2}\). Dividing this by \(t_{U}^{2}\) and assuming that the within-PSU sampling fraction, \(\bar{n} / N_{i}\), is negligible, we obtain the relvariance of \(\hat{t}_{pwr}\) as, approximately, \[\begin{equation} \label{eq:vpwr.2st} \small \frac{V\left(\hat{t}_{pwr} \right)}{t_{U}^{2} } \doteq \frac{B^{2} }{m} +\frac{W^{2} }{m\bar{n}} =\frac{\tilde{V}}{m\bar{n}} k \left[1+\delta \left(\bar{n}-1\right)\right] \end{equation}\]

with \(\small \tilde{V}=S_{U}^2/\bar{y}_{U}^2\), \(\small k=(B^2 + W^2)/\tilde{V}\), \[\begin{equation} \label{eq:B2.2st.pwr} \small B^{2} =\frac{S_{U1\left(pwr\right)}^{2} }{t_{U}^{2} } , \end{equation}\]

\[\begin{equation} \label{eq:W2.2st.pwr} \small W^{2} =\frac{1}{t_{U}^{2} } \sum _{i\in U}N_{i}^{2} \frac{S_{U2i}^{2} }{p_{i} }, \end{equation}\]

\[\begin{equation} \label{eq:delta.2st.pwr} \small \delta =B^{2} \left/ \left(B^{2} +W^{2} \right) \right. \end{equation}\]

Expression \(\small \eqref{eq:vpwr.2st}\) has the same form as \(\small \eqref{eq:deltavar.2st}\) but with different definitions of \(\small B^2\) and \(\small W^2\). Expression \(\small \eqref{eq:vpwr.2st}\) also has the interpretation of an srs variance of an unclustered variance, \(\small \tilde{V}/m\bar{n}\), times a design effect, \(\small k[1 + \delta(\bar{n}-1)]\), in the same way that \(\small \eqref{eq:deltavar.2st}\) did.

• Example. Between and within variance components in ppswr/srs design This example repeats the calculations in the example above for the variables in the Maryland area population. Assume that clusters will be selected proportional to the count of persons in each cluster. The function BW2stagePPS computes the population values of \(\small B^{2}\), \(\small W^{2}\), and \(\small \delta\) shown in \(\small \eqref{eq:B2.2st.pwr}\), \(\small \eqref{eq:W2.2st.pwr}\), and \(\small \eqref{eq:delta.2st.pwr}\) which are appropriate for ppswr sampling of clusters. The code for y1 using PSU or SSU as clusters is shown below. The variables, pp.PSU and pp.SSU, hold the one-draw probabilities \(\small p_{i}\) that appear in \(\small \eqref{eq:v2st}\):

 pp.PSU <- table(MDarea.popA$PSU) / nrow(MDarea.popA)
 pp.SSU <- table(MDarea.popA$SSU) / nrow(MDarea.popA)
 BW2stagePPS(MDarea.popA$y1, pp=pp.PSU, psuID=MDarea.popA$PSU)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#> 0.007878320 1.455574516 1.463117997 1.463452835 1.000228853 0.005383378
 BW2stagePPS(MDarea.popA$y1, pp=pp.SSU, psuID=MDarea.popA$SSU)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#>  0.03652265  1.42825477  1.46311800  1.46477742  1.00113417  0.02493392

The code for PSUs that are tracts and block groups is

 pp.trt <- table(MDarea.popA$TRACT) / nrow(MDarea.popA)
 pp.BG <- table(trtBG) / nrow(MDarea.popA)
 BW2stagePPS(MDarea.popA$y1, pp=pp.trt, psuID=MDarea.popA$TRACT)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#> 0.009458290 1.454058071 1.463117997 1.463516361 1.000272271 0.006462716
 BW2stagePPS(MDarea.popA$y1, pp=pp.BG, psuID=trtBG)
#>          B2          W2 unit relvar       B2+W2           k       delta 
#>  0.01643324  1.44797251  1.46311800  1.46440575  1.00088014  0.01122178

The between term when clusters are defined by PSU is about the same as when clusters are selected by srs because PSU’s all have the same size. With PSUs being either tracts or block groups in the ppswr/srswor design, the between term is much smaller than the within, compared to the results in the srs/srs example. For example, with y1 and srs sampling of tracts, \(\small B^2=0.2604\) but for pps sampling of tracts \(\small B^2=0.0091\).

When clusters are selected by srs, \(\small S_{U1}^{2}\) is the variance of the cluster totals around the average cluster total. In contrast, with pps sampling of clusters, \(\small S_{U1\left(pwr\right)}^{2}\) is the variance of the estimated population totals, \(\small t_{i} \left/ p_{i} \right.\) around the population total, \(\small t_{U}\). When clusters are selected with probability proportional to \(\small N_{i}\), then \(\small t_{i} \left/ p_{i} \right. = N_{i} \bar{y}_{Ui}\). If these one-cluster estimates of the population total are fairly accurate, as they are here, the \(\small B^{2}\) term can be quite small. This leads to much smaller values of \(\small \delta\) in pps sampling of clusters. This implies that the negative effect of clustering on the variance is lessened for a design that selects clusters with \(\small pp(N_i)\). This kind of comparison explains most practitioners’ preference for pps sampling of clusters, especially when the clusters vary in population size.