2.3 Theoretical foundations
As Heeringa et al. (2017) emphasize, the calculation of population totals and means, together with their variances, has been essential for the development of probability sampling theory and the proper interpretation of household survey results. Determining population totals is one of the fundamental pillars of survey analysis. Means, proportions, and ratios are all derived from these totals. A total is defined as the sum of a specific variable (for example, income or expenditure) across the entire population.
With household surveys, the analysis of numerical data frequently involves calculating descriptive statistics such as means, totals, and ratios, since these summarize the main characteristics of the population and serve as a basis for decision-making. These estimates can be calculated for the population as a whole or for specific subgroups, depending on the objectives of the research.
2.3.1 Point estimation
Once the objective of the sample design has been defined, the process of estimating the parameters of interest is carried out. For the purposes of this document, we begin with the estimation of total household income. For the estimation of totals in complex sample designs, strata will be denoted by the letter \(h\); primary sampling units, which are contained within strata, by the letter \(i\); while observed units (households or persons) will be denoted by the letter \(k\). Thus, the estimator of the total can be expressed as:
\[ \hat{t}_y = \sum_{h} \sum_{i} \sum_{k} w_{hik} \, y_{hik} \]
Here, \(w_{hik}\) corresponds to the expansion factor of unit \(k\), while \(y_{hik}\) corresponds to the observation of the variable of interest for that same unit. Calculating the estimate of the total and its associated variance can be complex, so it is necessary to consider different methodological approaches.
In terms of notation, and in order to simplify the writing throughout this document, unless otherwise indicated, \(\sum_{h}\) will denote the sum over all strata in the population; \(\sum_{i}\) will represent the sum over all primary sampling units selected in the sample within stratum \(h\); and \(\sum_{k}\) will correspond to the sum over all elements observed in the sample of UPM \(i\), belonging to stratum \(h\) of the population of interest.
2.3.2 Variance estimation
When working with household surveys, it is essential not only to obtain point estimates but also to quantify the uncertainty associated with those estimates. Understanding and estimating this uncertainty is an essential part of the analysis of household survey data. By applying appropriate methods, users can assess the precision of their estimates. There are several methods for estimating this precision and, with the support of modern software, these approaches can be implemented efficiently to support rigorous and reliable analyses. The main methods include:
Estimating equations: provide a flexible framework for estimating totals, means, ratios, and other parameters, as well as their corresponding variances, integrating a unified view of sampling theory (Binder, 1983).
Taylor linearization: consists of approximating complex nonlinear statistics through linear expressions and then estimating the variance of that approximation.
Ultimate cluster method: frequently used in surveys based on stratified multistage sampling. It is based on calculating variance from the differences between estimates obtained at the level of the primary sampling units (UPM). This method is often combined with Taylor linearization to estimate the variance of nonlinear statistics, such as means or ratios.
Bootstrap and other replication methods: are based on repeatedly taking subsamples from the observed data set, calculating estimates for each replicate, and then using the variability among these replicated estimates to infer the variance of the main estimator.
In the context of household surveys, many complex parameters can be formulated as linear combinations of other parameters \(f(\theta_1,\ldots,\theta_J)\). Specifically, let:
\[ f=f(\theta_1,\ldots,\theta_J) = \sum_{j=1}^J a_j \theta_j, \]
where the coefficients \(a_j\) are known constants and \(\theta_j\) represent population parameters. If \(\hat\theta_j\) is an estimator of \(\theta_j\), then a sample estimator of \(f\) is defined as:
\[ \hat{f} = \sum_{j=1}^J a_j \hat{\theta}_j, \]
Thus, its variance can be expressed as:
\[ Var(\hat{f}) = \sum_{j=1}^J a_j^2 Var(\hat{\theta}_j) \;+\; 2\sum_{j=1}^{J-1}\sum_{k>j}^J a_j a_k \, Cov(\hat{\theta}_j,\hat{\theta}_k). \]
This result is particularly important because it includes numerous cases of practical interest that will be developed throughout this document.
2.3.2.1 Estimating equations
Many population parameters can be expressed as solutions to estimating equations involving population totals. Although the technical details can be complex, the fundamental idea is that the same principles used to estimate totals can also be applied to variance estimation. This general framework makes the method simple and flexible, facilitating its implementation in specialized statistical software. A generic population estimating equation can be expressed as follows:
\[ \sum_{k\in U} z_k(\theta)=0, \]
where \(z_k(\cdot)\) is an estimating function evaluated for unit \(k\) and \(\theta\) represents the population parameter of interest. These equations provide a general framework for defining and calculating various population parameters, such as totals, means, and ratios. For example, for the population total, \(z_k(\theta)=y_k-\theta/N\) is defined, so that the estimating equation is \(\sum_{k\in U}(y_k-\theta/N)=0\). The solution to this equation leads to \(\theta=\sum_{k\in U} y_k = Y\), that is, to the population total.
Similarly, for the population mean, \(z_k(\theta)=y_k-\theta\) is used, and the equation \(\sum_{k\in U}(y_k-\theta)=0\) has as its solution \(\theta=\left(\sum_{k\in U} y_k\right)/N = \overline{Y}\), corresponding to the population mean. Likewise, in the case of ratios of totals, \(z_k(\theta)=y_k-\theta x_k\) is defined. Thus, the equation \(\sum_{k\in U}(y_k-\theta x_k)=0\) leads to the solution \(\theta=\dfrac{\sum_{k\in U} y_k}{\sum_{k\in U} x_k} = R\), which corresponds to the population ratio between the totals of variables \(y\) and \(x\). These formulations show how different parameters of interest can be expressed through estimating equations that share a common structure.
The idea of defining population parameters as solutions to estimating equations at the population level naturally leads to a general method for obtaining sample estimators. In this case, equations of the following form are used:
\[ \sum_{k\in s} w_k\, z_k(\theta)=0, \]
where \(w_k\) are the expansion factors and \(z_k(\theta)\) is the estimating function evaluated for each unit in the sample. Under probability sampling and assuming complete response, the sample sum \(\sum_{k\in s} d_k\, z_k(\theta)\) is unbiased with respect to its population analogue, which ensures that the solutions to these equations are consistent estimators of the population parameters.
2.3.2.2 Taylor linearization
This technique makes it possible to approximate the variance of nonlinear estimators. The procedure consists of applying a first-order Taylor expansion around the estimated parameter in order to replace the nonlinear estimator with a linear expression. This facilitates the calculation of variances in situations where exact formulas do not exist or their derivation is too complex. A consistent estimator of the variance, derived through Taylor linearization for solutions to sample estimating equations, can be expressed as:
\[ \widehat{Var}(\hat{\theta}) \approx [\hat{J}(\hat{\theta})]^{-1} \, \widehat{Var}_p \Bigg[\sum_{k\in s} w_k\, z_k(\hat{\theta})\Bigg] \, [\hat{J}(\hat{\theta})]^{-1} \]
where \(\hat{J}(\hat{\theta}) = \sum_{k\in s} w_k \left[ \frac{\partial z_k(\theta)}{\partial \theta} \right]_{\theta=\hat{\theta}}\). This result shows how Taylor linearization converts variance estimation for complex parameters into a problem of total estimation, which explains its widespread adoption in specialized software for survey analysis.
2.3.2.3 Ultimate cluster
The ultimate cluster method is a direct and robust approach for estimating the variance of totals in surveys that use stratified multistage cluster sampling designs. Proposed by Hansen et al. (1953), this method simplifies the complexity of multilevel designs by focusing only on the variation among primary sampling units (UPM). It is assumed that, within each sampling stratum, the UPM were selected independently with replacement (possibly with unequal probabilities), although in practice selection is usually carried out without replacement.
The method is based on the variation among statistics calculated at the UPM level. When applied correctly, it implicitly reflects any subsampling carried out within the UPM, allowing simpler but reliable variance estimates. It is especially useful in complex designs that include stratification and unequal selection probabilities for both UPM and lower-level units (households and individuals). The requirements for applying this method include the availability of unbiased estimates of totals for the variables of interest in each selected UPM. In addition, if the sample is stratified at the first stage, at least two UPM must be available per stratum, since this makes it possible to estimate adequately the variability among clusters within each stratum.
Thus, consider a multistage sampling design where \(n_h\) UPM are selected in stratum \(h\), (\(h=1,\dots,H\)). An estimate of the total of the variable of interest in UPM \(i\) of stratum \(h\) is:
\[ \hat{t}_{y_i} = \sum_{k\in s_{hi}} w_{hik} \ y_{hik} \]
Here, \(s_{hi}\) corresponds to the sample of units in UPM \(i\) of stratum \(h\); therefore, an unbiased estimator of the population total would be expressed as
\[ \hat{t}_{y} = \sum_{h}\sum_{i} \hat{t}_{y_i} \]
According to the ultimate cluster method, assuming that \(\hat{\bar t}_{y_h}=(1/{n_h}) \sum_{i} \hat{t}_{y_i}\) is the sample mean of the estimated UPM totals in stratum \(h\), an approximate estimator of the variance of \(\hat t_y\) is obtained through the following expression:
\[ \widehat{Var}(\hat{t}_y) = \sum_{h} \frac{n_h}{n_h-1} \sum_{i} (\hat{t}_{y_i} - \hat{\bar t}_{y_h})^2 \]
For more details, see Hansen, Hansen et al. (1953, p. 258) or Wolter & Wolter (2007). Although this method was originally proposed for calculating variances of total estimators, it can be combined with Taylor linearization or estimating equations to derive variances of other population parameters that can be formulated as solutions to estimating equations. This flexibility makes the method applicable to various contexts of household survey analysis.
A fundamental assumption of this technique is that, within each stratum, the UPM are selected independently and with replacement. In practice, most surveys select UPM without replacement, generating more efficient designs. Consequently, the variances calculated under the independence assumption are approximations of the true sampling variances. When the sampling fraction is small (for example, less than \(5\%\)), these approximations are usually sufficiently precise for use by national statistical offices or other analysts.
This method stands out for its simplicity and robustness, making it very attractive in practice. Although more sophisticated methods that consider all stages of the design may offer slightly more precise variance estimates, their application requires more detailed information and greater computational complexity. By contrast, the ultimate cluster method provides a reliable and efficient approximation, especially useful when estimating totals or means in household surveys. For a detailed analysis of the precision of this approximation and possible alternatives, see Särndal et al. (2003).
2.3.2.4 Bootstrap
In many cases, public survey microdata omit essential design information, such as identifiers for strata or primary sampling units (UPM), to protect respondent confidentiality. This omission limits users’ ability to calculate valid variances. In such situations, it is recommended that national statistical offices provide replication weights, which allow analysts to estimate standard errors correctly. Without these data, secondary users cannot reproduce the published standard errors or adequately account for the complex survey design.
Replication methods estimate variance by generating subsets of the original sample, calculating estimates for each one, and using the observed variability among these estimates to approximate the variance of the main estimator. They are particularly useful when information on strata or UPM is not available, a situation in which the ultimate cluster method cannot be applied.
For example, Bootstrap is a robust and versatile replication tool. Originally introduced by Efron (1979) for data that did not come from surveys, its most widely used adaptation for household surveys is the Rao-Wu-Yue Rescaling Bootstrap (Rao et al., 1992). This method fits optimally with stratified and multistage sampling designs and is widely used for variance estimation in complex surveys.
The procedure consists of generating many replicates of the original sample, simulating repeated draws from the population. Each replicate is constructed by creating additional columns of replication weights in the database, following this process:
- For each stratum, UPM are randomly selected with replacement; some may be repeated and others may not appear. Each selected UPM is included with all its observations. If the first-stage sample size in stratum \(h\) is greater than two (\(n_h > 2\)), the number of UPM selected per replicate is \(n_h - 1\).
- This process is repeated many times, usually hundreds, generating a large number of replicates. The number of times that a UPM \(i\) from stratum \(h\) appears in replicate \(r\) is denoted \(n_{hi}^{(r)}\), varying between 0 and \(n_h - 1\).
- From each replicate, new bootstrap weights are calculated for all units, reflecting how many times their UPM was selected. The weight of unit \(k\) in replicate \(r\) is calculated as:
\[ w_{hik}^{(r)} = w_{hik} \times \frac{n_h}{n_h - 1} \times n_{hi}^{(r)} \]
If the original weights include nonresponse or calibration adjustments, these must also be applied to each set of bootstrap weights.
When only Bootstrap replication weights are provided, analysts can estimate standard errors correctly, even without strata or UPM identifiers. For each replicate \(r\), the parameter of interest \(\hat{\theta}^{(r)}\) is calculated using the bootstrap weights \(w_{hik}^{(r)}\). The variance of the original estimator is approximated through the variability among all replicates:
\[ \widehat{Var}_B(\hat{\theta}) = \frac{1}{R} \sum_{r=1}^{R} \left(\hat{\theta}^{(r)} - \tilde{\theta}\right)^2, \quad \tilde{\theta} = \frac{1}{R} \sum_{r=1}^{R} \hat{\theta}^{(r)} \]
This approach ensures that the dispersion among replicates faithfully captures the uncertainty of the parameter. Bootstrap offers multiple advantages. Although it requires greater computational processing, it is effective for complex survey designs and makes it possible to estimate parameters that are difficult to calculate with traditional methods, such as medians or other nonlinear statistics. It is especially useful for analysts working with databases that lack strata and UPM identifiers but include replication weights.
The simplicity of the method facilitates its application even without specialized statistical software. However, most modern statistical packages already include procedures for applying Bootstrap and calculating variances, expanding its availability and robustness. Nevertheless, its use is not recommended in repeated surveys with overlapping samples or in situations with large sampling fractions and small sample sizes (Bruch et al., 2011).