2.1 Estimation of population means
In survey sampling, basic inferential procedures are developed for the estimation of finite population means. For the study variable \(y\), the population mean under stratification is given by
\[\begin{align} \mu_{y}= \frac{1}{N} \sum_{h=1}^H \sum_{i_h=1}^{N_h} y_{i_h} \,. \tag{2.1} \end{align}\]
Estimation of a disease prevalence is a special case of estimating a population proportion with a binary study variable \(y\). The design-unbiased Horvitz-Thompson estimator of \(\mu_y\) is given by
\[\begin{align}
\hat \mu_{y HT}= \frac{1}{N} \sum_{h=1}^H \sum_{i_h\in\mathcal{S}_h} w_{i_h} y_{i_h} \,.
\tag{2.2}
\end{align}\]
The stratum design weight \(w_{i_h} = \pi_{i_h}^{-1}\) is often interpreted as the number of units in the population represented by the unit \(i_h\) in the sample (Lohr 2010; Wu and Thompson 2020). The Horvitz-Thompson estimator for the population total \(T_y = \sum_{h=1}^H \sum_{i_h=1}^{N_h} y_{i_h}\) is given by \(\hat{T}_{y HT} = \sum_{h=1}^H \sum_{i_h\in\mathcal{S}_h} w_{i_h} y_{i_h}\), which is also called the expansion estimator. The design weight \(w_{i_h}\) is also called the inflation weight. The population size \(N\) is sometimes unknown to data users. An unbiased estimator of \(N\) is given by \(\hat{N} = \sum_{h=1}^H \sum_{i_h\in\mathcal{S}_h} w_{i_h}\). The resulting estimator of \(\mu_y\) is the so-called H\(\acute{\text{a}}\)jek estimator given by \(\hat{\mu}_{y H} = \hat{T}_{y HT} / \hat{N}\).
The theoretical design-based variance of the Horvitz-Thompson estimator given in (2.2) involves both the first-order and the second-order sample inclusion probabilities \(\pi_{i_h}\) and \(\pi_{i_h i_h'}\). Under stratified sampling, the stratum samples \(\mathcal{S}_h\), \(h=1,\cdots,H\) are independent. The general theoretical variance formula is given by \[\begin{align*} \mbox{Var}(\hat{\mu}_{y HT}) &= N^{-2} \sum^H_{h=1} \mbox{Var}\left( \sum_{i_h\in\mathcal{S}_h} \frac{y_{i_h}}{\pi_{i_h}} \right) \\ & = N^{-2} \sum^H_{h=1} \left[ \sum_{i_h=1}^{N_h} \sum_{i_h^\prime = 1}^{N_h} \Big(\pi_{i_h i_h^\prime }- \pi_{i_h} \pi_{i_h^\prime } \Big) \frac{y_{i_h}}{\pi_{i_h}} \frac{y_{i_h}}{\pi_{i_h}} \right] \,. \end{align*}\] The conventional unbiased variance estimator for the Horvitz-Thompson estimator is given by \[\begin{align} \widehat{\mbox{Var}}(\hat{\mu}_{y HT}) &= N^{-2} \sum^H_{h=1}\left[ \sum_{i_h\in\mathcal{S}_h} \sum_{i_h'\in\mathcal{S}_h} \frac{\pi_{i_h i_h^\prime }- \pi_{i_h} \pi_{i_h^\prime }}{\pi_{i_h i_h' }} \frac{y_{i_h}}{\pi_{i_h}} \frac{y_{i_h}}{\pi_{i_h}} \right] \,. \tag{2.3} \end{align}\]
In practice, complex survey datasets, such as the CLSA datasets used in this paper, usually do not provide the joint inclusion probabilities \(\pi_{i_h i_h' }\) which are required for computing the variance estimator given in (2.3). Most statistical software packages for survey data analyses use approximate variance estimators to bypass this difficulty.
When sampled units are drawn with-replacement, with selection probabilities \(z_{i_h}\), \(i_h = 1,\cdots, N_h\) for each selection, the Hansen-Hurwitz estimator (Hansen and Hurwitz 1943) of \(\mu_y\) has the same algebraic form of the Horvitz-Thompson estimator if we let \(\pi_{i_h} = n_h z_{i_h}\). The unbiased variance estimator for the Hansen-Hurwitz estimator is given by \[\begin{align} \hat{V}_0 = N^{-2} \sum^H_{h=1} \frac{1}{n_h(n_h-1)} \sum_{i_h \in\mathcal{S}_h} \left(n_h w_{i_h}y_{i_h} - \sum_{i_h \in\mathcal{S}_h} w_{i_h}y_{i_h} \right)^2 \; . \tag{2.4} \end{align}\] The variance estimator \(\hat{V}_0\) does not involve second-order inclusion probabilities and provides a good approximation to the variance estimator given in ((2.3)) if the sampling fractions \(f_h = n_h/N_h\) are small for the original without-replacement survey design. When the sampling fractions are not small, an ad hoc adjustment to (2.4) is to apply the finite population correction factor \(1-f_h\) within each stratum. The resulting variance estimator is given by \[\begin{align} \hat{V} = N^{-2} \sum^H_{h=1} \frac{1-f_h}{n_h(n_h-1)} \sum_{i_h \in\mathcal{S}_h} \left(n_h w_{i_h}y_{i_h} - \sum_{i_h \in\mathcal{S}_h} w_{i_h}y_{i_h} \right)^2 \; . \tag{2.5} \end{align}\]
The variance estimator \(\hat{V}\) given in (2.5) is exactly design-unbiased for stratified simple random sampling. For general stratified unequal probability sampling, the performance of \(\hat{V}\) varies depending on the original survey design. There exist other approximate variance formulas not involving second-order inclusion probabilities and performing better for certain designs. See, for instance, (Haziza, Mecatti, and Rao 2008) for further details. The variance estimator \(\hat{V}\) is the default option for most survey packages, including \(\texttt{R,\ SAS,\ SPSS}\) and \(\texttt{Stata}\).
Stratified multi-stage sampling can use approximate variance estimators similar to \(\hat{V}\) if sampling fractions for the first stage clusters are small within each stratum. For cases where \(N\) is unknown and the H\(\acute{\text{a}}\)jek estimator \(\hat{\mu}_{y H}\) is used, the variance estimator \(\hat{V}\) given in (2.5) needs to be modified with \(N\) being replaced by \(\hat{N}\) and the study variable \(y_{i_h}\) being substituted by the residual variable \(e_{i_h} = y_{i_h} - \hat{\mu}_{y H}\) for computing the variance estimator. Further details can be found in (Wu and Thompson 2020).