B.1 Two Possible Inferences
In household surveys, it is useful to distinguish three objects. First, the finite population, denoted by \(U = \{1,2,\ldots,N\}\), which contains all units of interest. Second, the study variable \(y_k\), observed or defined for each unit \(k \in U\). Third, the sample design \(p(s)\), which assigns probabilities to the possible samples \(s \subset U\). In the notation used throughout this document, a unit can be identified as \(k\) within cluster or PSU \(i\), belonging to stratum \(h\), so weighted estimators are often written as sums of the form
\[ \hat{t}_y = \sum_h \sum_i \sum_k w_{hik} \ y_{hik} \]
where \(w_{hik}\) represents the expansion factor or sampling weight associated with the observed unit. In sampling theory, the values \(y_k\) are considered fixed once the population has been defined. Randomness is not in the values of the variable, but in the inclusion indicators
\[ I_k = \begin{cases} 1, & \text{if } k \in s,\\ 0, & \text{if } k \notin s, \end{cases} \]
which induce the inclusion probability \(\pi_k = \Pr_p(I_k=1)\) and, in turn, the corresponding sampling weight, defined as \(w_k = \frac{1}{\pi_k}\).
An estimator that ignores the sampling weights may not adequately represent the target population, especially when the sample comes from a complex design. In design-based inference, the population is considered fixed and the source of randomness comes from the sample selection mechanism; therefore, its central elements are the inclusion probabilities \(\pi_k\) and the sampling weights \(w_k\).
By contrast, model-based inference considers the values of the variable of interest across the population as realizations of a stochastic process, usually denoted by \(\xi\). Thus, while model-based inference studies properties such as \(E_{\xi}(Y_k)\) and \(Var_{\xi}(Y_k)\), design-based inference is oriented toward incorporating the characteristics of the sample plan into the estimators, so that they are unbiased \(E_{p}(\hat{t}_y) = t_y\).
This distinction makes it possible to understand why, even when there is a reasonable model for the variable of interest, the sample design must be incorporated explicitly if the goal is to obtain valid inferences for the target population. To understand the distinction between model-based inference and design-based inference, we start with a simple example, adapted from Binder (2011). Suppose that \(N=100\) independent realizations of a Bernoulli variable with parameter \(\theta=0.3\) are generated, where \(\theta\) represents the expected proportion of unemployed people in a population generated by a model. If \(Y_k \sim Bernoulli(\theta)\), then the expectation under the model is \(E_{\xi}(Y_k) = \theta\), while the variance under the model is \(Var_{\xi}(Y_k) = \theta(1-\theta)\).
In the specialized literature, the model \(\xi\) is called a superpopulation model, because it generates finite populations. Thus, for a finite population generated by this model, the population mean is \(\bar{Y}_U = \frac{1}{N}\sum_{k \in U}Y_k\). Note that, under the model \(\xi\), this mean is unbiased for \(\theta\), since
\[ E_{\xi}(\bar{Y}_U) = E_{\xi}\left(\frac{1}{N}\sum_{k \in U}Y_k\right) = \frac{1}{N}\sum_{k \in U}E_{\xi}(Y_k) = \theta. \]
The following Monte Carlo simulation reproduces this process. In each repetition, a complete population of size \(N=100\) is generated and its mean is calculated. The average of these population means should approximate \(\theta\).
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_model <- 1000
model_estimates <- rep(NA, n_sim_model)
for (sim in seq_len(n_sim_model)) {
outcome <- rbinom(population_size, 1, theta)
model_estimates[sim] <- mean(outcome)
}
cbind(
theta,
expected_model_mean = mean(model_estimates),
bias = mean(model_estimates) - theta
)## theta expected_model_mean bias
## [1,] 0.3 0.302 0.00225The code sets a seed so that the exercise is reproducible. It then generates n_sim_model independent populations using rbinom(). The output reports the true parameter theta, the average of the simulated means expected_model_mean, and the Monte Carlo bias. As expected under the IID Bernoulli model, expected_model_mean is very close to 0.3 and the bias is close to zero. The small differences occur because 1000 simulations are used, not an infinite number of repetitions.
In sampling theory, by contrast, the characteristics of interest are fixed population parameters. If a person is unemployed, their status is considered a fixed value of the finite population. What is random is the mechanism through which the sample is selected. The population total is defined as \(t_y = \sum_{k \in U} y_k\), while the population mean is \(\bar{y}_U = \frac{t_y}{N}\). If a sample \(s\) is observed, the Horvitz-Thompson estimator of the total is
\[ \hat{t}_{y} = \sum_{k \in s}\frac{y_k}{\pi_k} = \sum_{k \in s}w_k y_k. \]
As mentioned in previous chapters, when the population size \(N\) is unknown, the population mean can be estimated using the weighted ratio estimator, defined as
\[ \hat{\bar{y}} = \frac{\sum_{k \in s}w_k y_k} {\sum_{k \in s}w_k} \]
These estimators incorporate the inclusion probabilities into their functional form, and this is the central difference from the simple unweighted average of a sample. To illustrate the impact of the design, suppose that the previous population is divided into \(N_I\) PSUs. PSU \(i\) has size \(N_i\), total \(t_{y_i}=\sum_{k \in U_i}y_{ik}\), and mean \(\bar{y}_i=t_{y_i}/N_i\). If \(n_I\) PSUs are selected with probabilities proportional to size, then their inclusion probabilities take the following form:
\[ \pi_{Ii} = \frac{n_I N_i}{N} \]
When all members of the selected PSUs are observed, the weighted estimator of the population mean, which is approximately unbiased, can be written as
\[ \hat{\bar{y}} = \frac{1}{N}\sum_{i \in S_I}\frac{t_{y_i}}{\pi_{Ii}} = \frac{1}{N}\sum_{i \in S_I}\frac{N_i\bar{y}_i}{n_I N_i/N} = \frac{1}{n_I}\sum_{i \in S_I}\bar{y}_i. \]
This expression illustrates why, when the sampling design is complex (in this case, probability proportional to size selection), the average of the cluster means provides unbiasedness. By contrast, the simple average of the people observed in the selected PSUs would be:
\[ \bar{y}_s = \frac{\sum_{i \in S_I}t_{y_i}} {\sum_{i \in S_I}N_i}, \]
which treats the sample as if it were self-weighting. This estimator is not unbiased when the selection probabilities are unequal or when cluster size is related to the variable of interest. The following simulation fixes a finite population and repeatedly selects households with probabilities proportional to size using S.piPS() from the TeachingSampling package (Gutiérrez, 2020). In each sample, two estimators are calculated: design_estimates, which corresponds to \(\hat{\bar{y}}\), and simple_estimates, which corresponds to the simple average of the people observed in the selected PSUs, \(\bar{y}_s\). In addition, these estimators are compared with the model parameter theta_population, corresponding to \(\theta = 0.3\).
library(TeachingSampling)
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_sampling <- 1000
design_estimates <- simple_estimates <- rep(NA, n_sim_sampling)
cluster_size <- rep(2:6, each = 5)
n_clusters <- length(cluster_size)
cluster_id <- rep(seq_len(n_clusters), cluster_size)
size_center <- weighted.mean(cluster_size, cluster_size)
prob_person <- theta + rep((cluster_size - size_center) / 12, cluster_size)
outcome <- rbinom(population_size, 1, prob_person)
theta_population <- mean(outcome)
cluster_totals <- tapply(outcome, cluster_id, sum)
cluster_means <- tapply(outcome, cluster_id, mean)
sampled_cluster_count <- floor(n_clusters * 0.3)
for (sim in seq_len(n_sim_sampling)) {
pps_sample <- S.piPS(sampled_cluster_count, cluster_size)
sampled_clusters <- pps_sample[, 1]
sampled_cluster_size <- cluster_size[sampled_clusters]
sampled_cluster_total <- cluster_totals[sampled_clusters]
sampled_cluster_mean <- cluster_means[sampled_clusters]
design_estimates[sim] <- mean(sampled_cluster_mean)
simple_estimates[sim] <-
sum(sampled_cluster_total) / sum(sampled_cluster_size)
}
cbind(
theta_population,
expected_pps = mean(design_estimates),
bias_pps = mean(design_estimates) - theta_population,
expected_simple = mean(simple_estimates),
bias_simple = mean(simple_estimates) - theta_population
)## theta_population expected_pps bias_pps expected_simple bias_simple
## [1,] 0.3 0.304 0.00368 0.358 0.0585The first block of code generates households of unequal sizes, between 2 and 6 people. To make the design effect visible, the individual probability of unemployment varies smoothly with household size, but is centered so that the expected population average remains close to theta. Then theta_population is calculated, which is the actual proportion of unemployment in that specific population. In each repetition, S.piPS() selects PSUs with probabilities proportional to size; then the mean of the means of the selected PSUs and the simple average of the observed people are calculated. The output compares both estimators with theta_population. The Monte Carlo bias of expected_pps is close to zero, while that of expected_simple is not, showing the result of ignoring the design in this exercise.