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.