Appendix B Finite Population Inference
A fundamental question in survey analysis is to identify which probability measure governs the inference. In classical statistical inference, the available observations are usually assumed to be independent and identically distributed (IID) random variables. Under this approach, randomness comes from the model that generates the data. However, as Kish (1987) notes, this assumption does not adequately describe information from surveys with complex sample designs, where observations are selected through stratification, clustering, and unequal inclusion probabilities.
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.
B.2 The Integration of Both Inferences
Nevertheless, in many applications these two inferential frameworks are not presented as mutually exclusive approaches, but as complementary components of the same problem. In particular, the sample may be considered to have been obtained through a complex probability design while, at the same time, the values of the variable of interest respond to a superpopulation model. This formulation leads to a double-inference framework, in which uncertainty does not come from a single source, but from the combination of the mechanism that generates the population values (superpopulation model) and the mechanism that determines which units are observed (sampling design).
Formally, two probability measures are involved: \(\xi\), associated with the process that generates the values \(Y_k\), and \(p(\cdot)\), associated with the sample design that selects the sample \(s\). Therefore, the properties of the estimators must be evaluated while simultaneously considering both sources of variation. Under this framework, the combined expectation can be expressed as
\[ E_{\xi p}(\hat{\theta}) = E_{\xi}\{E_p(\hat{\theta}\mid U)\}, \]
where the behavior of the estimator is first evaluated under the design, conditional on the finite population generated, and then averaged over the superpopulation model. In this context, the full joint likelihood may be intractable, especially when the inclusion probabilities depend on variables related to the response. The solution proposed by Pfeffermann (1993) is the maximum pseudo-likelihood (MPV) technique, which weights each unit’s contribution by the inverse of its inclusion probability.
B.2.1 Maximum Likelihood Method
The maximum likelihood method starts from a random sample \(y_1,\ldots,y_n\) generated by a distribution with density or probability function \(f(y;\theta)\). If the observations are independent and identically distributed (IID), the likelihood function is defined as
\[ L(\theta) = \prod_{k=1}^{N} f(y_k;\theta). \]
Because products can be difficult to manipulate, the log-likelihood is used:
\[ \ell(\theta) = \sum_{k=1}^{N}\log f(y_k;\theta) \]
The maximum likelihood estimator \(\hat{\theta}\) is the value of \(\theta\) that maximizes \(\ell(\theta)\). If the function is differentiable, this value satisfies the score equations:
\[ U(\theta) = \frac{\partial \ \ell(\theta)}{\partial \ \theta} = \sum_{k=1}^{N} u_k(\theta) = 0 \]
where \(u_k(\theta) = \frac{\partial}{\partial\theta}\log f(y_k;\theta)\) is the contribution of unit \(k\) to the total score. For a Bernoulli distribution with parameter \(\theta\), the probability function is \(f(y_k;\theta) = \theta^{y_k}(1-\theta)^{1-y_k}\), for \(y_k \in \{0,1\}\). The log-likelihood is
\[ \ell(\theta) = \sum_{k=1}^{N} \left[ y_k\log(\theta) +(1-y_k)\log(1-\theta) \right] \]
By differentiating, setting equal to zero, and solving, the MV estimator is obtained, defined by:
\[ \hat{\theta}_{MV} = \frac{1}{N}\sum_{k=1}^{N}y_k \]
Thus, the sample proportion is the natural estimator of the Bernoulli parameter. The key point is that this result depends on the assumption that the observations have the same distribution and do not come from a design with unequal probabilities.
On the other hand, in a multiple linear regression model with normal errors, we have \(\mathbf{y} \sim N(\mathbf{X}\boldsymbol{\beta},\sigma^2 I)\). In this case, the log-likelihood, omitting constants, can be written as
\[ \ell(\boldsymbol{\beta},\sigma^2) = -\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})'(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \]
Maximizing this expression with respect to \(\beta\) is equivalent to minimizing the sum of squared residuals. Therefore, the maximum likelihood estimator of \(\beta\) coincides with the ordinary least squares estimator \(\hat{\boldsymbol{\beta}}_{MV} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\). This result is important because it shows that many common statistical procedures can be interpreted as maximum likelihood under specific model assumptions. However, if the data come from a complex survey, the matrix \(\mathbf{X}'\mathbf{X}\) does not by itself reflect the sample selection mechanism.
B.2.2 Maximum Pseudo-Likelihood Method
When the sample design is complex, the observed data do not necessarily satisfy the assumptions of independence and identical distribution. Some units may represent many people in the population, and others only a few. In addition, two units selected within the same cluster may be correlated. In this context, directly using the ordinary log-likelihood can produce estimates that describe the sample, but not the target population.
Maximum pseudo-likelihood, developed from the results of Binder (1983), incorporates the sampling weights into the log-likelihood. For observed units \(k \in s\), it is defined as
\[ \ell_p(\theta) = \sum_{k \in s} w_k \log f(y_k;\theta) \]
The weight \(w_{hik}=1/\pi_{hik}\) indicates how many population units the observation \(y_{hik}\) represents. Therefore, the contribution of a unit with a low inclusion probability is amplified in the pseudo-likelihood. The resulting estimating equations are
\[ U_p(\theta) = \sum_{k \in s} w_k \ u_k(\theta) = 0 \]
MPV estimation preserves the logic of maximum likelihood, but replaces the ordinary score with a weighted score. The solution \(\hat{\theta}_{MPV}\) is the value that makes the weighted sum of individual contributions equal to zero. For a Bernoulli distribution, the pseudo-log-likelihood is
\[ \ell_p(\theta) = \sum_{k \in s} w_k \left[ y_k\log(\theta) +(1-y_k)\log(1-\theta) \right]. \]
By differentiating and setting equal to zero, we obtain:
\[ \frac{\partial l_p(\theta)}{\partial \theta} = \sum_{k \in s} w_k \left[ \frac{y_k}{\theta} -\frac{1-y_k}{1-\theta} \right] =0 \]
The solution is:
\[ \hat{\theta}_{MPV} = \frac{\sum_{k \in s}w_k y_k} {\sum_{k \in s}w_k} = \hat{p}_d \]
Therefore, for a binary variable, MPV leads to the weighted estimator of the proportion, equivalent to the Hájek estimator. This result directly connects pseudo-likelihood theory with the proportion estimators presented in the chapters on categorical variables.
Likewise, for a multiple linear regression model, the weighted pseudo-log-likelihood implies minimizing a weighted sum of squares:
\[ (\mathbf{y}-\mathbf{X}\beta)'\mathbf{W}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}), \]
where
\[ \mathbf{W}= \begin{pmatrix} w_1 & 0 & \cdots & 0\\ 0 & w_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & w_n \end{pmatrix}. \]
The solution is the weighted least squares estimator:
\[ \hat{\boldsymbol{\beta}}_{MPV} = (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}\mathbf{X}'\mathbf{W}\mathbf{y} \]
This result naturally extends the ordinary estimator to the context of complex surveys: each observation contributes to the model fit in proportion to its sampling weight. Nevertheless, for inference it is not enough to weight the point estimate; the variance must also be estimated correctly, taking strata, clusters, and weights into account.
The following nested Monte Carlo simulation illustrates double inference. In the outer loop, a new population is generated under the Bernoulli model. In the inner loop, many samples are selected by design and the estimator that incorporates the design is compared with the simple sample average.
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_population <- 100
expected_design_estimates <- expected_simple_estimates <- population_means <-
rep(NA, n_sim_population)
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)
sampled_cluster_count <- floor(n_clusters * 0.3)
for (population_sim in seq_len(n_sim_population)) {
outcome <- rbinom(population_size, 1, prob_person)
population_means[population_sim] <- mean(outcome)
cluster_totals <- tapply(outcome, cluster_id, sum)
cluster_means <- tapply(outcome, cluster_id, mean)
design_estimates <- simple_estimates <- rep(NA, 100)
for (sample_sim in seq_len(100)) {
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[sample_sim] <- mean(sampled_cluster_mean)
simple_estimates[sample_sim] <-
sum(sampled_cluster_total) / sum(sampled_cluster_size)
}
expected_design_estimates[population_sim] <- mean(design_estimates)
expected_simple_estimates[population_sim] <- mean(simple_estimates)
}
cbind(
theta,
mean_population = mean(population_means),
expected_pps = mean(expected_design_estimates),
bias_pps = mean(expected_design_estimates) - theta,
expected_simple = mean(expected_simple_estimates),
bias_simple = mean(expected_simple_estimates) - theta
)## theta mean_population expected_pps bias_pps expected_simple bias_simple
## [1,] 0.3 0.304 0.302 0.00248 0.332 0.0317The code reproduces the two sources of randomness. The outer loop represents the model \(\xi\), because it generates new finite populations with rbinom(). The inner loop represents the design \(p\), because multiple samples are selected for each population with S.piPS(). The output includes mean_population, which is the average of the generated population means; expected_pps, which summarizes the behavior of the estimator that incorporates the design; and expected_simple, which summarizes the simple average of the selected sample. When the design is incorporated correctly, the Monte Carlo bias with respect to theta should be small. The comparison with expected_simple shows why survey-based inference requires weights and design variances, even when starting from a simple statistical model.
In summary, ordinary maximum likelihood is appropriate when observations can be treated as IID under a probabilistic model. In complex surveys, maximum pseudo-likelihood offers a natural extension: it preserves the structure of the estimating equations, but weights each individual contribution by its sampling weight. Thus, statistical models are fitted in a way that is more coherent with the population the survey seeks to represent.