The Generalized Multivariate Difference Estimator (GMDe) in Sample Survey Statistics: Formulation, Properties, and Comparative Analysis
I. Introduction to the Generalized Multivariate Difference Estimator (GMDe)
A. Conceptual Overview and Rationale
The Generalized Multivariate Difference estimator (GMDe) has emerged as a significant development in sample survey statistics, particularly for complex and longitudinal surveys. Its primary objective is to enhance the precision of population estimates for multiple study variables by leveraging information from various auxiliary variables.1 These auxiliary variables are typically correlated with the study variables and are often less costly or more readily available to measure, thereby offering a pathway to reduce uncertainty in the final estimates.3 In contemporary survey environments, such as national forest inventories 2 or extensive longitudinal studies, researchers often have access to a rich array of auxiliary data. This can include information from remote sensing technologies, administrative records, outputs from deterministic models, or data from previous survey waves or different sample phases. The GMDe provides a systematic framework for integrating these diverse sources of information to produce more statistically efficient estimates for the key variables under investigation.
The capacity of GMDe to synthesize information from disparate sources positions it as a methodological advancement towards data fusion within the domain of survey statistics. It allows for the coherent integration of estimates derived from various systems or measurement technologies, such as combining field plot re-measurements with remotely sensed data for disturbance detection and time-series predictions from ecosystem models.1 The recursive version, rGMDe, further facilitates the integration of time-series of multivariate population estimates from multiple systems.4 This capability extends beyond the simple use of auxiliary variables; it implies a mechanism for combining different estimates or data layers, each with its own potential imperfections and error structures. Such data fusion is critical for constructing a more comprehensive and potentially more accurate understanding of population parameters, especially in dynamic environments or when tracking changes over time. This makes GMDe a potentially powerful tool for statistical agencies and research programs that manage and analyze multiple, evolving data streams concerning the same or related populations.
B. Origins and Relationship to Classical Difference Estimators
The GMDe is not an entirely novel concept in its foundational principles but rather a multivariate extension and generalization of the classical univariate difference estimator, a technique well-documented by Hansen, Hurwitz, and Madow in 1953.1 The classical difference estimator is a cornerstone in traditional survey sampling theory, designed to improve the estimation of a population mean or total by utilizing an auxiliary variable for which the population mean or total is known.5 The work of Hansen, Hurwitz, and Madow is considered foundational in the development of sample survey methods and theory.7 GMDe takes this established concept and expands its applicability to scenarios involving multiple study variables and, crucially, multiple auxiliary variables simultaneously.
The characterization of the classical difference estimator family as "long-neglected" in some of the literature pertaining to GMDe 10 suggests that its core principle—using differences relative to auxiliary information to reduce the variance of an estimator—has held untapped potential, particularly in multivariate contexts. The efficiency of the classical univariate difference estimator hinges on the strength of the correlation between the study variable and the auxiliary variable. A strong positive correlation means that the difference between the study variable and a scaled auxiliary variable will have a smaller variance than the study variable itself. If this fundamental principle is sound and effective in the univariate case, its extension to a multivariate framework, where multiple such correlations across various study and auxiliary variables can be exploited concurrently, should theoretically lead to substantial gains in statistical efficiency. The historical "neglect" of a full multivariate generalization might have been due to the computational or theoretical complexities that arise when moving from a single study and auxiliary variable to multiple dimensions. Modern computational capabilities and the theoretical framework provided by GMDe aim to overcome these past limitations, thereby revitalizing and significantly expanding a proven statistical idea for contemporary data-rich environments.
II. Mathematical Formulation and Core Principles of GMDe
A. Notation and Setup: Study Variables and Auxiliary Residuals
The Generalized Multivariate Difference estimator (GMDe) operates on a vector of initial estimates for M study variables and J auxiliary residuals. Let y be the M×1 vector of study variables of primary interest to the analyst, and let x represent information related to J auxiliary variables. The GMDe typically starts with a set of design-based estimates for population parameters, such as population totals or means. For instance, one might begin with Horvitz-Thompson estimates derived from a probability sample.
Let T^y​ be an M×1 vector of initial design-based estimates for the population totals of the M study variables.
Let T^x∗​ be a J×1 vector of initial design-based estimates related to the J auxiliary variables.
A crucial aspect of GMDe is its use of "auxiliary residuals" rather than necessarily the raw auxiliary variables or their known population totals. An auxiliary residual, denoted e^x​, is a J×1 vector representing deviations or differences. For example, an auxiliary residual could be the difference between a population estimate for an auxiliary variable from a Phase-One sample and a corresponding estimate from a Phase-Two sample, or the difference between a sample estimate of an auxiliary variable and a known population total if available, or residuals from a model predicting the auxiliary variables.1
The GMDe process commences with an (M+J)×1 vector, T^, of these initial, typically design-consistent, estimates:
T^=(T^y​e^x​​)
Associated with T^ is its (M+J)×(M+J) design-based covariance matrix, V(T^), which is partitioned conformably:
$$ \mathbf{V}(\hat{\mathbf{T}}) = \begin{pmatrix} \mathbf{V}{yy} & \mathbf{V}{yx} \ \mathbf{V}{xy} & \mathbf{V}{xx} \end{pmatrix} $$
Here, Vyy​ is the M×M covariance matrix of T^y​, Vxx​ is the J×J covariance matrix of e^x​, and Vyx​ (with Vxy​=Vyx′​) is the M×J covariance matrix between T^y​ and e^x​.1 The concept of auxiliary residuals allows GMDe to be flexible; it can operate even when true population totals for auxiliary variables are unknown, relying instead on differences between various estimates of these auxiliary quantities.
B. The GMDe Transformation and Coefficient Matrix
The GMDe is fundamentally a simple linear transformation applied to the initial vector of estimates. The adjusted estimates for the study variables, T^y,GMDe​, are obtained as:
T^y,GMDe​=T^y​−Be^x​
where B is an M×J matrix of coefficients (or weights).1 This matrix B embodies the "model" or the empirically determined relationship used to adjust the initial estimates of the study variables using the information contained in the auxiliary residuals. The product Be^x​ represents an M×1 vector of adjustments.
More generally, the GMDe can be expressed as a linear transformation of the larger vector T^ using an M×(M+J) coefficient matrix, let's call it C. If C is partitioned as $$, where IM​ is the M×M identity matrix, then T^y,GMDe​=CT^. The constants in this coefficient matrix are computed from sample survey estimates, specifically from the covariance matrix V(T^).1
C. The Minimum Variance Criterion
The core principle for determining the "optimal" coefficient matrix B is the minimization of the variance of the resulting GMDe estimates, T^y,GMDe​, for each of the M study variables.1 The variance of T^y,GMDe​ is given by:
$$ \mathbf{V}(\hat{\mathbf{T}}{y,GMDe}) = \mathbf{V}(\hat{\mathbf{T}}y - \mathbf{B} \hat{\mathbf{e}}x) = \mathbf{V}{yy} - \mathbf{B} \mathbf{V}{xy} - \mathbf{V}{yx} \mathbf{B}' + \mathbf{B} \mathbf{V}_{xx} \mathbf{B}' $$
To find the matrix B that minimizes the variance for each study variable (i.e., minimizes the diagonal elements of V(T^y,GMDe​)), standard multivariate theory leads to the solution:
Bopt​=Vyx​Vxx−1​
This optimal M×J coefficient matrix, Bopt​, is analogous to the matrix of regression coefficients in a multivariate linear regression of T^y​ on e^x​. It is parameterized using the design-based estimates of the covariance matrix partitions Vyx​ and Vxx​ (or, equivalently, the correlation matrix between study variables and auxiliary residuals, and the covariance/correlation matrix of the auxiliary residuals themselves).10
The degree of variance reduction achieved by GMDe depends directly on the strength of the (multiple) correlation between the study variables and the auxiliary residuals, as captured by these covariance matrices.10 If the auxiliary residuals are highly correlated with the estimation errors in T^y​, then Bopt​e^x​ can effectively predict and subtract out a substantial portion of this error, leading to a T^y,GMDe​ with lower variance.
The "optimality" of GMDe, in this context, is conditional on the sample-estimated covariance matrix V(T^). This implies a critical dependency: the quality and accuracy of the initial design-based estimates of Vyx​ and Vxx​ are paramount for the effective performance of GMDe.10 If these initial covariance (and thus correlation) estimates are imprecise or biased, the resulting "optimal" coefficient matrix Bopt​ will be suboptimal, potentially leading to diminished variance reduction or, in worst-case scenarios, even variance inflation compared to the initial estimates T^y​. This underscores the importance of a well-designed initial survey phase that can yield reliable estimates of these second-order statistics.
III. The Recursive GMDe (rGMDe): Addressing Computational Challenges and Enhancing Robustness
A. Motivation for the Recursive Approach
A significant practical challenge arises in the direct application of the GMDe formulation, specifically in the computation of the optimal coefficient matrix Bopt​=Vyx​Vxx−1​. This computation requires the inversion of Vxx​, the J×J covariance matrix of the auxiliary residuals.2 In many contemporary survey applications, the number of potential auxiliary variables (J) can be very large. With a large J, or when auxiliary variables exhibit strong multicollinearity, the matrix Vxx​ can become rank-deficient (singular) or ill-conditioned (near-singular). In such cases, the inverse Vxx−1​ either does not exist or is numerically unstable, rendering the calculation of Bopt​ infeasible or unreliable.2 This computational hurdle is a major limitation for the standard GMDe, especially given the increasing availability and use of high-dimensional auxiliary datasets (e.g., from remote sensing or administrative records) where multicollinearity is common.
B. The rGMDe Algorithm: Recursive Scalar Inversions
To overcome the infeasibility of inverting a potentially large and rank-deficient Vxx​ matrix, the recursive GMDe (rGMDe) algorithm was developed.2 The central innovation of rGMDe is that it replaces the direct inversion of the J×J matrix Vxx​ with a sequential process involving J scalar inversions, one for each of the J auxiliary residuals, processed one at a time.2 This recursive approach ensures that the rGMDe is always computationally feasible, irrespective of the dimensionality (J), rank, or condition number of the full auxiliary residual covariance matrix Vxx​.10 This makes the GMDe framework practically applicable across a much broader spectrum of survey scenarios, particularly those characterized by a large number of auxiliary variables.
C. Handling Rank-Deficient and Ill-Conditioned Covariance Matrices
The recursive nature of rGMDe, by sidestepping the need for direct matrix inversion of Vxx​, inherently manages situations where this matrix is rank-deficient or ill-conditioned.2 If, during a particular recursion step, an auxiliary residual is found to have zero or near-zero variance (after accounting for previously processed residuals through orthogonalization), this indicates a linear dependency or severe multicollinearity. In such instances, that specific auxiliary residual is deemed infeasible for that step and is simply omitted from the adjustment process for that recursion, preventing numerical instability.2 This adaptive characteristic ensures the numerical stability and reliability of the rGMDe algorithm, which is crucial for an estimator intended for complex, high-dimensional datasets where multicollinearity among auxiliary variables is a frequent occurrence. The ability of rGMDe to bypass the problematic matrix inversion is a direct cause of its enhanced robustness in these high-dimensional settings. This makes rGMDe particularly well-suited for integrating information from "big data" sources, where the number of potential auxiliary variables (J) can be exceptionally large, and inter-correlations are common.3
D. Stepwise Selection of Auxiliary Residuals
The rGMDe algorithm incorporates an efficient stepwise selection procedure within each of its J recursions.2 At each step j (from 1 to J), the algorithm evaluates the remaining J−j+1 auxiliary residuals (that have not yet been fully processed) to identify the one that is most effective in reducing the variance of the population estimates for the study variables. The selection can be guided by pre-defined heuristic weights provided by the analyst to prioritize variance reduction for specific, more critical study variables.2 This stepwise approach aims to construct a sensible and efficient weight matrix by iteratively incorporating the most informative auxiliary residual at each stage. This feature adds an element of intelligent variable selection, focusing computational effort on the auxiliary information that provides the greatest marginal benefit in terms of variance reduction. It also helps to mitigate the risk of incorporating spurious correlations by prioritizing stronger, more impactful relationships 3, and can proceed without requiring intervention from an experienced statistician during the computation.3
E. Recursive Orthogonalization
A key component of the rGMDe algorithm is the recursive orthogonalization of the auxiliary residuals.2 As each auxiliary residual is selected and processed in a given recursion, its effect is, in a sense, removed from the remaining auxiliary residuals for subsequent recursions. This is achieved by computing orthogonal weights such that the auxiliary residual being processed at the current step becomes numerically independent of (uncorrelated with) those residuals processed in prior steps and those yet to be processed in future steps, within the context of their contribution to adjusting the study variables. This orthogonalization process is critical for preventing the "double counting" of information that would occur if correlated auxiliary variables were used naively. It ensures that each selected auxiliary residual contributes uniquely to the variance reduction, leading to more stable and interpretable weights, especially when the original auxiliary variables exhibit multicollinearity. This combination of stepwise selection and recursive orthogonalization in rGMDe functions somewhat analogously to a forward stepwise regression procedure. However, instead of optimizing for model fit (e.g., R2), the rGMDe optimizes for variance reduction of the survey estimates, which is a primary concern in survey sampling. It sequentially extracts the independent contribution of each selected auxiliary variable towards this goal.
F. Inequality Constraints: Min-Max Criteria, Outlier Mitigation, and Coherent Covariance Matrix
A significant enhancement in the rGMDe framework is its ability to relax the strict minimum variance criterion to incorporate various inequality constraints.2 These constraints are applied during each recursion to ensure the practical utility and statistical soundness of the final estimates. They include:
1.	Minimum-Maximum Criteria: The resulting population estimates for study variables can be bounded by analyst-defined minimum and maximum plausible values.2 This prevents estimates from falling into unrealistic ranges.
2.	Outlier Mitigation: The algorithm includes mechanisms to identify and mitigate the influence of suspiciously large or influential auxiliary residuals.2 This is achieved by modifying the minimum variance coefficients to down-weight the impact of potential outliers, enhancing the robustness of the estimator.
3.	Coherent Covariance Matrix: Constraints are imposed to ensure that the final estimated M×M covariance matrix for the M study variables is "coherent".2 This means ensuring that all estimated variances are non-negative and all estimated correlations fall within the logical range of ±1. While this might involve a slight sacrifice in the theoretical minimum variance, it ensures that the resulting covariance matrix is statistically valid and interpretable.3
4.	Equality Constraints: The rGMDe can also impose equality constraints, for example, to ensure "additivity" in statistical tables (e.g., estimates for sub-domains sum to the estimate for the parent domain) or consistency across different levels of a geospatial hierarchy.10 Such constraints can improve the overall accuracy and consistency of the set of estimates, much like calibration estimators achieve.
These constraints collectively enhance the robustness and practical applicability of the rGMDe estimates, ensuring they are not only efficient but also sensible and internally consistent.
IV. Key Statistical Properties of GMDe and rGMDe
A. Variance of GMDe Estimates and Factors Influencing Efficiency
The central aim of the GMDe is the reduction of variance in population estimates for the study variables, thereby improving their statistical efficiency.1 The extent of this variance reduction is directly linked to the strength of the correlation between the study variables and the auxiliary residuals used in the adjustment.1 Specifically, if T^y,GMDe​=T^y​−Bopt​e^x​ with Bopt​=Vyx​Vxx−1​, the variance of the GMDe estimator is:
$$ \mathbf{V}(\hat{\mathbf{T}}{y,GMDe}) = \mathbf{V}{yy} - \mathbf{V}{yx} \mathbf{V}{xx}^{-1} \mathbf{V}_{xy} $$
This can be rewritten in terms of the multiple correlation matrix. Let Ryx​ be the M×J matrix of correlations between the M study variables and J auxiliary residuals, and Rxx​ be the J×J correlation matrix of the auxiliary residuals. The term Vyx​Vxx−1​Vxy​ represents the reduction in variance. The higher the squared multiple correlation between each study variable and the set of auxiliary residuals, the greater the variance reduction. Thus, the selection of highly correlated auxiliary variables is crucial for maximizing the efficiency gains from GMDe.
B. Bias Considerations, particularly in the Estimated Covariance Matrix
While the GMDe point estimator, T^y,GMDe​, aims to be design-consistent (and potentially unbiased under certain conditions if the initial estimates are unbiased and B is considered fixed), a significant issue arises with the estimation of its variance-covariance matrix, V(T^y,GMDe​). The GMDe methodology, in its standard formulation, tends to produce negatively biased estimates of this M×M covariance matrix.2 This means that the estimated variances are, on average, smaller than the true variances, and confidence intervals constructed using these underestimated variances will be too narrow, leading to an overstatement of precision and potentially incorrect statistical inferences.
This underestimation of variance is a common problem in survey statistics when weights or adjustment factors are themselves estimated from the sample data but are subsequently treated as fixed constants in the variance calculation. In the case of GMDe, the coefficient matrix Bopt​ is a function of Vyx​ and Vxx​, which are estimated from the sample. Thus, Bopt​ is a random matrix. If the variance formula for T^y,GMDe​ does not account for the variability arising from the estimation of Bopt​, it will lead to an underestimation of the true variance.2 The GMDe estimates are conditional on the specific sample realized, as a different sample would yield different estimates for Vyx​ and Vxx​, and consequently a different Bopt​ matrix.10 To address this negative bias in variance estimation, resampling methods such as the bootstrap are recommended.2 The bootstrap, by re-estimating Bopt​ in each bootstrap sample, can more accurately capture the total variability.
C. Underlying Assumptions and Dependence on Design-Based Estimates
The GMDe procedure commences with design-based estimates of population parameters (e.g., totals or means for study variables and auxiliary residuals) and, critically, their design-based variance-covariance matrix.10 The Horvitz-Thompson estimator is often cited as an example of the initial design-based estimator. The performance and accuracy of the GMDe are heavily reliant on the quality of these initial inputs. In particular, the accuracy of the estimated M×J population correlation (or covariance) matrix between the study variables and the auxiliary residuals, derived from this design-based component, is crucial for the parameterization of the "optimal" coefficient matrix B.10 If these initial correlation estimates are poor or inaccurate, the resulting GMDe adjustments will be suboptimal, potentially negating the benefits of the method.
The nature of GMDe concerning model dependency warrants careful consideration. It has been described variously as "model-free" 1 and "model-based".10 A more precise characterization appears to be that of a hybrid or model-assisted approach where the "assisting model" is empirically derived. The process begins with design-based estimates, which are model-free in the sense that their validity relies on the probability design rather than an assumed superpopulation model relating the study variables to auxiliary variables. However, GMDe then introduces a model-based component: the linear adjustment Be^x​. The parameters of this adjustment model (the elements of B) are not pre-specified based on an assumed superpopulation relationship but are estimated empirically from the design-based covariance structure.10 Thus, GMDe uses the data to build its own optimal linear adjustment model, conditional on the observed relationships in the sample. This distinguishes it from purely model-based estimators that assume a superpopulation model from the outset, and also from purely design-based estimators that make no use of such modeling for adjustment. It is "model-assisted" in that a model is used to improve efficiency, but the model itself is data-driven from design-valid inputs.
V. Comparative Analysis: GMDe/rGMDe versus Alternative Estimators
The GMDe positions itself as a pioneering alternative to several conventional estimators in survey statistics.2 A comparative analysis with these alternatives, particularly other multivariate estimators, is essential to understand its specific contributions and potential advantages.
A. Generalized Regression Estimator (GREG)
The Generalized Regression Estimator (GREG) is a widely used model-assisted technique in survey sampling.13
1. Univariate vs. Multivariate Nature:
GMDe is inherently a multivariate estimator, designed to simultaneously handle multiple study variables (M>1) and produce estimates of their population parameters along with their full M×M covariance matrix.2 This is a key distinction from the standard GREG, which is primarily a univariate estimator, typically applied to one study variable at a time.2 While GREG can be applied repeatedly to each of several study variables, the estimation of the covariance matrix for the resulting vector of GREG estimates is not an intrinsic part of the standard GREG procedure and often requires separate, more complex considerations, such as using weighted residuals within a suitable variance estimation framework for the specific sampling design.13 Some research has explored multivariate extensions or alternatives to GREG, such as multivariate ratio estimators 17, further highlighting that the standard GREG is not directly multivariate in the same way as GMDe. A comparative study noted that GMDe "allows several study variables to be estimated at the same time," contrasting with GREG's typical application.18
2. Handling of Auxiliary Information and Model Assumptions:
GMDe utilizes population estimates for auxiliary variables (often in the form of auxiliary residuals, e^x​) and their estimated variance-covariance matrix Vxx​.1 Crucially, GMDe does not require that the population totals (or means) of the auxiliary variables be known constants without error; it can operate with auxiliary information that is itself subject to sampling or prediction error.2
In contrast, the traditional GREG estimator assumes that the population totals of the auxiliary variables (the X-totals) are known precisely, for example, from a census or administrative register.13 GREG is model-assisted; it employs a working linear model (superpopulation model) that specifies the relationship between the study variable y and the auxiliary variables x to improve the efficiency of the Horvitz-Thompson estimator. The GREG estimator is designed to be approximately design-unbiased and design-consistent even if this working model is misspecified, although its efficiency benefits are greatest when the model is a good approximation of reality.13 The model primarily serves to determine the optimal form of the adjustment to the Horvitz-Thompson estimator.
3. Optimization Goals and Efficiency:
The optimization goal of GMDe is to determine the coefficient matrix B that minimizes the variance of the GMDe estimator for each of the M study variables, based on the sample-estimated covariances (or correlations) between the study variables and the auxiliary residuals.1
The GREG estimator's primary "calibration" property is to ensure that the sample-weighted sum of the auxiliary variables equals their known population totals. This calibration property, when the auxiliary variables are correlated with the study variable, leads to a reduction in the variance of the estimator for the study variable's total or mean.6 From a modeling perspective, GREG seeks to find the best coefficient estimates for the assumed linear model between the survey data and the auxiliary data, and then uses these to adjust the Horvitz-Thompson estimator.18
4. High-Dimensional Covariates and Two-Level GREG:
GREG estimators have been extended to handle situations with high-dimensional auxiliary covariates (p is large). This may involve techniques such as LASSO (Least Absolute Shrinkage and Selection Operator) for variable selection and coefficient estimation, or principal component analysis (PCA) to reduce the dimensionality of the auxiliary information before applying GREG.13 Specialized GREG estimators, such as two-level GREGs, have also been developed to ensure consistency between estimates at different hierarchical levels (e.g., person and household) in complex surveys.23
The rGMDe, by its design, is well-suited for high-dimensional auxiliary information (large J). Its recursive algorithm, which involves scalar inversions and stepwise selection of auxiliary residuals, allows it to accommodate hundreds or even thousands of auxiliary variables without the computational instability associated with inverting very large covariance matrices.2 This inherent capability for high dimensionality might offer a more direct and computationally stable approach compared to GREG in scenarios with a vast number of auxiliary variables. While GREG can be adapted, it might require pre-processing steps like PCA or regularization, which introduce their own set of assumptions and complexities. rGMDe's integrated stepwise selection and orthogonalization are designed to navigate high-dimensional spaces directly as part of its estimation procedure.
B. Mixed Estimators (ME) and Model-Based Approaches
Mixed Estimators (ME) represent a class of estimators that often combine information from different sources or statistical models, frequently incorporating random effects to account for unexplained heterogeneity or hierarchical data structures.
1. Conceptual Similarities and Differences:
GMDe is described as a type of composite estimator.2 It is also noted to be closely related to the concept of an empirical Best Linear Unbiased Predictor (BLUP) and the Kalman filter update equations, which are common in time series analysis and state-space models.1 While GMDe starts with design-based estimates, its adjustment phase is guided by an empirically derived model, leading to its classification as a "model-based method" in some contexts.10
Mixed Estimators are prominently used in Small Area Estimation (SAE), where the goal is to produce reliable estimates for domains with small sample sizes by "borrowing strength" from related areas or auxiliary data.25 These often involve specifying models with random area effects. The Fay-Herriot (FH) model is a well-known area-level mixed model that links direct survey estimates for small areas to area-specific covariates via a linear model with random area effects.29 Unit-level mixed models, such as the Battese-Harter-Fuller (BHF) model, relate unit-level survey responses to unit-level covariates and include random area effects.33 Both FH and BHF models are used to produce Empirical Best Linear Unbiased Predictors (EBLUPs) of small area means or totals.
2. Computational Aspects (Matrix Inversion vs. Recursion):
A key practical difference highlighted is that traditional Mixed Estimators often require the inversion of covariance matrices, which can become computationally intensive or unstable if these matrices are large or ill-conditioned (e.g., when dealing with many small areas or many parameters).2 The rGMDe algorithm, with its recursive sequence of scalar inversions, is specifically designed to avoid this limitation, making it more scalable for problems with a large number of auxiliary variables (which, in some GMDe applications, might play a role analogous to many small effects or areas).2
3. Relationship to BLUP, Fay-Herriot, and other Unit/Area-Level Models (e.g., BHF):
The characterization of GMDe as an empirical BLUP 2 is significant. BLUP, in general, provides the best (minimum mean squared error) linear unbiased predictor of random effects or linear combinations involving random effects in a linear mixed model framework.37 Estimators derived from FH and BHF models are EBLUPs, where variance components are estimated from the data.
The connection suggests that GMDe shares foundational statistical principles with these established mixed models. However, GMDe's operational approach differs. Instead of explicitly modeling random effects and estimating variance components as in FH or BHF models, GMDe directly adjusts initial design-based estimates using an empirically optimized linear transformation based on estimated covariances. This might imply that GMDe offers a more direct, and potentially less assumption-laden, route to "borrowing strength" or improving estimates when the primary goal is the estimation of general population parameters (which could include aggregates over small areas) rather than the explicit prediction of area-specific random effects themselves. GMDe's framework might be more flexible if the specific distributional or structural assumptions of particular mixed models (like FH or BHF) do not hold well, although it might be less efficient if those model assumptions are indeed correct and the model is well-specified.
C. Calibration Estimators
Calibration estimators adjust the original design weights of sample units so that the weighted sums of auxiliary variables in the sample conform to known population totals for these auxiliary variables.42
1. Weight Adjustment Principles:
Calibration fundamentally involves modifying sampling weights. A distance measure (e.g., chi-square) between the initial design weights (dk​=1/πk​) and the final calibrated weights (wk​) is minimized subject to a set of calibration equations, typically of the form ∑k∈s​wk​xk​=XTotal​, where XTotal​ are known population totals for a vector of auxiliary variables xk​.6
GMDe, in its core formulation, adjusts population estimates rather than directly modifying individual sample unit weights.1 However, the rGMDe variant can impose equality constraints.10 For instance, it can constrain the sum of small-area GMDe estimates to match a pre-determined estimate for a larger area. This functional outcome—achieving consistency with known totals or superordinate estimates—is very similar to the aim of calibration estimators, although the mechanism (adjusting estimates via a coefficient matrix vs. adjusting unit weights) differs.
2. Relationship to GREG and Asymptotic Equivalence:
The GREG estimator is a prominent member of the family of calibration estimators; in fact, GREG can be derived through the calibration framework using a specific distance function.42 Many estimators derived using different reasonable distance functions within the calibration framework are asymptotically equivalent to the GREG estimator, meaning they share similar large-sample properties, including variance.43
3. Model-Calibration (Wu & Sitter) and Extensions:
The model-calibration approach, introduced by Wu and Sitter (2001), extends classical calibration by allowing the use of more complex working models (e.g., nonlinear, nonparametric, or even machine learning models like neural networks) to predict the study variable yk​ using auxiliary variables xk​.45 The fitted values, y^​k​, from such a model are then used in the calibration equations, typically ∑k∈s​wk​y^​k​=∑k∈U​y^​k​ (if population-level fitted values can be obtained) or calibrating to the sum of y^​k​ over the sample if only sample-level covariates are used to fit the model that generates y^​k​. This allows calibration to leverage more sophisticated relationships than the simple linear ones implicit in standard GREG.
GMDe's adjustment is based on auxiliary residuals (e^x​) rather than fitted values of the study variables (y^​k​) in the manner of model-calibration. While both use models, the nature of the model and how it enters the estimation differs.
4. GMDe's Constraint Capabilities in Relation to Calibration:
The capacity of rGMDe to impose equality constraints 10—for example, to ensure that estimates for sub-populations sum to a known or estimated total for the entire population, or that small area estimates are consistent with estimates for larger geographical aggregations—makes it functionally analogous to calibration estimators in achieving such consistencies. This suggests a convergence in the desired properties of modern survey estimators: the ability to produce a coherent set of estimates that align with external information or internal additivity requirements. While calibration achieves this by modifying weights and GMDe achieves it by constraining the final adjusted estimates, both methodologies offer paths to this important objective, albeit from different starting points and with different underlying mechanics.
D. Other Relevant Multivariate Estimators (Brief Overview)
To fully situate GMDe, it is useful to acknowledge other multivariate estimation techniques used in or relevant to survey statistics:
●	Multivariate Regression (mvreg): This is a standard statistical technique for modeling the relationship between multiple dependent variables and one or more predictor variables simultaneously.48 It typically assumes multivariate normality of residuals and is used for predictive modeling rather than direct population total estimation in a design-based survey context, though its coefficients might inform model-assisted estimators.48
●	Multivariate Random-Effects Meta-Analysis: Used to synthesize results from multiple studies when each study reports multiple, correlated outcomes. It accounts for both within-study and between-study correlations to produce pooled estimates for each outcome and their covariance.50 This is relevant when survey estimates from different cycles or related surveys need to be combined.
●	Multivariate Logistic-Assisted Estimators: These are specialized for categorical or multinomial response variables in complex surveys, using logistic regression models in a model-assisted framework to estimate totals or proportions for multiple categories.52
●	Principal Component Based Estimators: In situations with a large number of auxiliary variables (high-dimensional x), PCA can be used as a dimension reduction technique. The principal components, which are linear combinations of the original auxiliary variables, can then be used in regression estimation or calibration procedures to avoid multicollinearity and improve stability.13
GMDe's approach of directly adjusting a vector of design-based estimates using an empirically optimized transformation based on auxiliary residuals is distinct from these methods, though it shares the goal of leveraging multivariate information.
The following table provides a comparative overview of GMDe/rGMDe against GREG, Mixed Estimators, and Calibration Estimators.
Table 1: Comparative Overview of Key Estimators
Feature	GMDe/rGMDe	GREG (Standard)	Mixed Estimators (e.g., Fay-Herriot for SAE)	Calibration Estimators (General)	
Primary Nature	Multivariate	Univariate	Univariate or Multivariate (model-dependent)	Typically Univariate (weights applied to single y) but can ensure consistency for multiple X's	
Auxiliary Info Requirement	Population estimates of auxiliary parameters & their covariance; residuals often used.	Known population totals for auxiliary variables (X-totals).	Unit or area-level auxiliary data; explicit model structure (e.g., random effects).	Known population totals for auxiliary variables (X-totals).	
Core Mechanism	Linear transformation of initial estimates using an empirically optimized coefficient matrix.	Adjustment of Horvitz-Thompson estimator based on a working linear model.	Prediction based on a mixed-effects model (EBLUP).	Adjustment of design weights to match known X-totals.	
Handling Multiple Study Vars	Simultaneous estimation with full covariance matrix.	Applied one-by-one; covariance requires separate consideration.	Can be modeled jointly if multivariate mixed model is specified.	Weights are common, but each y-total estimated separately.	
Key Strengths	Multivariate covariance; robust to rank-deficiency (rGMDe); flexible with auxiliary data types.	Design-consistency; model-robustness (for point estimate); widely understood.	Efficient for small areas if model holds; explicit modeling of heterogeneity.	Enforces consistency with known totals; can improve precision.	
Key Limitations/Assumptions	Potential negative bias in variance estimation; relies on quality of initial covariance estimates.	Assumes known X-totals; less direct for multivariate covariance.	Relies on model specification; design properties can be complex (e.g., BHF).	Can produce extreme weights; choice of distance function.	
Robustness Features (rGMDe/General)	rGMDe: outlier mitigation, inequality/equality constraints, handles multicollinearity, censors spurious correlations.	Asymptotically design-unbiased even if working model is misspecified.	Depends on model assumptions; some ME are robust to certain deviations.	Some methods can control for extreme weights.	
VI. Advantages, Disadvantages, and Application Contexts of GMDe/rGMDe
A. Strengths
The GMDe, particularly in its recursive (rGMDe) formulation, offers several compelling advantages for survey estimation:
1.	Multivariate Capability: It is inherently designed to simultaneously estimate parameters for multiple study variables and their full variance-covariance matrix, providing a comprehensive view of their interrelationships.1
2.	Robustness of rGMDe: The recursive algorithm exhibits significant robustness. It effectively handles rank-deficient or ill-conditioned auxiliary covariance matrices by avoiding direct matrix inversion.2 It incorporates mechanisms to mitigate the influence of outliers in auxiliary data and can censor spurious or weak correlations to prevent them from unduly influencing the estimates.2 Furthermore, the imposition of inequality and equality constraints ensures that the final estimates are coherent (e.g., non-negative variances, valid correlations, additivity).2
3.	Flexibility with Auxiliary Data: GMDe is not restricted to using known population totals for auxiliary variables. It can effectively utilize estimates of auxiliary parameters, including residuals from other models or differences between estimates from various sources (e.g., multi-phase samples).1 This allows for the integration of diverse data types, such as remote sensing data, outputs from deterministic models, administrative records, and "big data" sources.2
4.	High-Dimensional Data Handling: The rGMDe algorithm is capable of accommodating a very large number of auxiliary variables—potentially hundreds or thousands—due to its recursive scalar inversion process and stepwise variable selection.2
5.	Adaptability: The framework is adaptable to evolving survey objectives and measurement technologies within long-term survey programs, allowing for the integration of new data sources over time.10
6.	Conceptual Simplicity: Despite its power, the underlying principle of GMDe is a linear transformation that extends the classical, well-understood difference estimator, making its core concept relatively straightforward.1
B. Limitations
Despite its strengths, GMDe also has limitations that users must consider:
1.	Potential Bias in Covariance Estimation: A primary concern is the tendency for GMDe to produce negatively biased estimates of the M×M variance-covariance matrix of the study variable estimates.2 This underestimation of variance can lead to overly optimistic precision measures and invalid confidence intervals. Resampling techniques like the bootstrap are recommended to mitigate this bias.2
2.	Nascent Nature of the Method: GMDe is a relatively recent development in survey statistics. While its theoretical properties are being established, the authors themselves note that some assertions require further rigorous challenge and more extensive empirical comparisons with conventional methods in diverse and complex survey settings.11
3.	Dependence on Accurate Initial Correlation Estimates: The performance of GMDe, particularly the degree of variance reduction achieved, is highly dependent on the accuracy of the initial design-based estimates of the covariance (and thus correlation) matrix between the study variables and the auxiliary residuals.10 Inaccurate or unstable initial correlation estimates will lead to a suboptimal adjustment model and may reduce or even negate the estimator's benefits.
C. Suitability for Complex Surveys, Longitudinal Data, and Integration of Modern Data Sources
GMDe is explicitly designed for application in complex sample surveys, including those with longitudinal components where data are collected over time.1 Its ability to process large vectors of population estimates, which can represent different time periods, various domains (including small areas), and cells within multifaceted statistical tables, underscores its suitability for large-scale institutional surveys like National Forest Inventories.2
The framework's capacity to integrate auxiliary information from a wide array of modern data sources—such as remote sensors, deterministic model outputs, administrative records, and "big data"—makes it particularly relevant in the current data-rich environment.3 This positions GMDe, especially the rGMDe variant, as a sophisticated "meta-estimator" or an "integrating estimator." It can synthesize information from various prior estimation processes or diverse data streams within a unified system. By starting with "prior estimates" and their associated covariances 1 and being capable of integrating "time-series of multivariate population estimates from multiple systems" 4, GMDe offers a powerful tool for statistical agencies and large research programs. These organizations often manage and seek to reconcile multiple data sources to produce the most coherent and precise understanding of population characteristics. The ability to handle auxiliary residuals, rather than requiring precisely known auxiliary population totals, further enhances its applicability in real-world scenarios where perfect auxiliary information is rare.
The following table summarizes key characteristics of GMDe and rGMDe.
Table 2: Summary of GMDe and rGMDe Characteristics
Feature	Description / Mathematical Essence / Implication	
Basic Principle	Multivariate extension of classical difference estimator; improves study variable estimates using correlated auxiliary variable information.	
Mathematical Form	T^y,GMDe​=T^y​−Be^x​; a linear transformation of initial estimates.	
Variance Objective	Minimizes variance for each study variable estimate, conditional on sample-estimated correlations/covariances.	
rGMDe - Matrix Inversion	Replaces inversion of J×J auxiliary covariance matrix (Vxx​) with J recursive scalar inversions.	
rGMDe - Rank-Deficiency Handling	Inherently handles rank-deficient/ill-conditioned Vxx​ by avoiding its direct inversion; infeasible residuals in recursion are omitted.	
rGMDe - Stepwise Selection	Selects the most effective auxiliary residual at each recursion step to maximize variance reduction.	
rGMDe - Orthogonalization	Recursively orthogonalizes auxiliary residuals to prevent double-counting of information from correlated auxiliaries.	
rGMDe - Inequality Constraints	Imposes min-max bounds on estimates, mitigates outlier influence, ensures coherent final covariance matrix (e.g., non-negative variances).	
rGMDe - Equality Constraints	Can enforce additivity (e.g., small area estimates sum to larger area totals), similar to calibration.	
Variance Estimation Issue	Standard GMDe variance formula tends to underestimate true variance (negative bias) because B is estimated. Bootstrap recommended for bias reduction.	
Key Assumptions / Dependencies	Starts with design-based estimates and their covariance matrix. Accuracy heavily depends on the quality of initial design-based correlation estimates between study and auxiliary variables.	
VII. Conclusion and Future Directions
A. Summary of GMDe/rGMDe's Position in Survey Statistics
The Generalized Multivariate Difference Estimator (GMDe), and particularly its recursive variant (rGMDe), represents a notable advancement in the field of sample survey statistics. It builds upon the foundational principles of the classical difference estimator, extending it to a multivariate framework capable of simultaneously improving estimates for multiple study variables by leveraging information from numerous, potentially diverse, auxiliary sources. GMDe can be best characterized as a model-assisted estimator where the "assisting model" (the coefficient matrix for adjustment) is empirically derived from design-based estimates of covariances and correlations, rather than being pre-specified from superpopulation assumptions.
The development of rGMDe has been crucial in addressing the practical computational challenges associated with inverting large and potentially ill-conditioned covariance matrices of auxiliary variables, a common issue in high-dimensional data settings. Through its recursive scalar inversions, stepwise selection of auxiliary residuals, orthogonalization procedures, and the incorporation of inequality and equality constraints, rGMDe offers a robust, flexible, and computationally feasible approach. Its key advantages lie in its inherent multivariate nature, its capacity to handle high-dimensional and complex auxiliary information (including estimates rather than only known totals), and its features for outlier mitigation and ensuring coherent outputs. These characteristics make GMDe/rGMDe particularly well-suited for complex, longitudinal surveys and for integrating data from modern sources like remote sensing, administrative records, and model predictions.
B. Potential Areas for Further Research and Validation
Despite its promising features, GMDe is a relatively nascent methodology, and further research is warranted to fully establish its properties and broaden its applicability. As noted by its proponents, the assertions regarding GMDe's performance, especially in comparison to well-established conventional methods, require continued rigorous challenge and validation across a wider array of complex survey designs and diverse data scenarios.11
Specific areas for future investigation include:
1.	Extensive Empirical Comparisons: More comprehensive simulation studies and real-world applications are needed to compare the performance (bias, variance, mean squared error) of GMDe/rGMDe against GREG, various Mixed Estimators (including unit-level and area-level SAE models), and different Calibration Estimators. These comparisons should cover a range of conditions, such as varying sample sizes, different levels of correlation between study and auxiliary variables, the presence of outliers, and complex sampling designs (e.g., multi-stage, informative).
2.	Variance Estimation Properties: While the bootstrap method is recommended for addressing the negative bias in the GMDe variance estimator, further theoretical and empirical research into the properties of this bootstrap variance estimator specifically for GMDe would be beneficial. This includes assessing its bias, stability, and coverage properties for confidence intervals under different conditions. Alternative analytical variance estimators that better account for the variability of the estimated coefficient matrix could also be explored.
3.	Diagnostic Tools: The development of diagnostic tools tailored for GMDe/rGMDe applications would aid practitioners in assessing the quality of the initial correlation estimates, the influence of individual auxiliary variables, the impact of constraints, and the overall goodness-of-fit of the GMDe adjustment.
4.	Sensitivity Analysis: Systematic studies on the sensitivity of GMDe to inaccuracies in the initial design-based covariance estimates and to the choice of parameters for constraint application (e.g., min-max bounds, outlier detection thresholds) would be valuable.
5.	Broader Application Domains: While GMDe has seen application in areas like forestry 2, exploring its utility in other fields that rely on complex survey data (e.g., official economic statistics, health surveys, social surveys) could demonstrate its versatility and highlight new areas for methodological refinement.
6.	Software Implementation: Wider availability of user-friendly and well-documented software to implement GMDe/rGMDe would facilitate its adoption and further testing by the broader survey statistics community. The mention of an R-function 2 and potential SQL implementation 11 are steps in this direction.
In conclusion, the GMDe framework, especially with the enhancements offered by rGMDe, provides a powerful and adaptable tool for survey statisticians seeking to improve the precision and coherence of multivariate population estimates by effectively integrating diverse sources of auxiliary information. Continued research and validation will further solidify its place within the array of modern survey estimation techniques.
Works cited
1.	(PDF) Generalized Multivariate Difference Estimator (GMDe) - ResearchGate, accessed May 16, 2025, https://www.researchgate.net/publication/353411938_Generalized_Multivariate_Difference_Estimator_GMDe
2.	Generalized Multivariate Difference Estimator (GMDe): The Recursive Algorithm (rGMDe) - ResearchGate, accessed May 16, 2025, https://www.researchgate.net/profile/Raymond-Czaplewski/publication/376983859_Generalized_Multivariate_Difference_Estimator_GMDe_The_Recursive_Algorithm_rGMDe/links/659ac3fb0bb2c7472b39fcbb/Generalized-Multivariate-Difference-Estimator-GMDe-The-Recursive-Algorithm-rGMDe.pdf
3.	Generalized Multivariate Difference Estimator (GMDe): The Recursive Algorithm (rGMDe), accessed May 16, 2025, https://www.researchgate.net/publication/376983859_Generalized_Multivariate_Difference_Estimator_GMDe_The_Recursive_Algorithm_rGMDe
4.	Generalized Multivariate Difference Estimator (GMDe) (1 book series) Hardcover Edition, accessed May 16, 2025, https://www.amazon.com/dp/B0CVG8CVN3?binding=hardcover&ref=dbs_dp_awt_sb_pc_thcv
5.	1987: VARIANCE ESTIMATION USING AUXILIARY INFORMATION UNDER A ONE UNIT PER STRATUM SAMPLE DESIGN - Proceedings of the Survey Research Methods Section, accessed May 16, 2025, http://www.asasrms.org/Proceedings/papers/1987_131.pdf
6.	Connections between survey calibration estimators and semiparametric models for incomplete data - PubMed Central, accessed May 16, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC3699889/
7.	The 2009 Morris Hansen Lecture: The Care, Feeding, and Training of Survey Statisticians - Washington Statistical Society, accessed May 16, 2025, https://washstat.org/hansen/2009_JOS_Lohr.pdf
8.	12. Sample Survey Methods and Theory, accessed May 16, 2025, https://academic.oup.com/jrsssa/article-pdf/158/2/350/49759892/jrsssa_158_2_350.pdf
9.	1997: OPTIMAL CLUSTERING IN MULTI-STAGE SAMPLES - Proceedings of the Survey Research Methods Section, accessed May 16, 2025, http://www.asasrms.org/Proceedings/papers/1997_133.pdf
10.	Generalized Multivariate Difference Estimator (GMDe): The Recursive Algorithm (rGMDe), accessed May 16, 2025, https://www.amazon.com/Generalized-Multivariate-Difference-Estimator-GMDe/dp/B0CVGDMS77
11.	(PDF) Generalized Multivariate Difference Estimator (GMDe): The ..., accessed May 16, 2025, https://www.researchgate.net/publication/383221694_Generalized_Multivariate_Difference_Estimator_GMDe_The_Recursive_Algorithm_rGMDe
12.	GMDe | Environmetrika, accessed May 16, 2025, https://environmetrika.org/gmde
13.	Generalized Regression Estimators with High-Dimensional Covariates - PMC, accessed May 16, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC7313320/
14.	Review of how the generalized regression estimators contribute to estimating the financial and economic data with missing observations under unequal probability sampling | Emerald Insight, accessed May 16, 2025, https://www.emerald.com/insight/content/doi/10.1108/ajeb-12-2023-0134/full/html
15.	Replication Methods for Two-phase Sampling - CRAN, accessed May 16, 2025, https://cran.r-project.org/web/packages/svrep/vignettes/two-phase-sampling.html
16.	www.asasrms.org, accessed May 16, 2025, http://www.asasrms.org/Proceedings/papers/2000_016.pdf
17.	On the Improvement of Multivariate Ratio Method of Estimation in Sample Surveys by Calibration Weightings - SciSpace, accessed May 16, 2025, https://scispace.com/pdf/on-the-improvement-of-multivariate-ratio-method-of-1eon6cis47.pdf
18.	Comparison Study of the Generalized Multivariate Difference ..., accessed May 16, 2025, https://www.causeweb.org/usproc/eusr/2022/virtual-posters/8
19.	A Short Note on Extended GREG Family Estimators for Domains - Official Statistics, accessed May 16, 2025, https://officialstatistics.wordpress.com/wp-content/uploads/2010/05/bok10.pdf
20.	A Comparison of Model-Assisted Estimators, With and Without Data-Driven Transformations of Auxiliary Variables, With Application to Forest Inventory - Frontiers, accessed May 16, 2025, https://www.frontiersin.org/journals/forests-and-global-change/articles/10.3389/ffgc.2021.764495/full
21.	Calibration Based on Principal Components - DiVA portal, accessed May 16, 2025, https://www.diva-portal.org/smash/get/diva2:574097/FULLTEXT02
22.	Principal Component Regression with Survey Data. Application on the French Media Audience 1 Introduction 2 Model-assisted and ca - ISI, accessed May 16, 2025, https://isi-web.org/sites/default/files/import/files-2011/950593.pdf
23.	(PDF) A Two‐level GREG Estimator for Consistent Estimation in ..., accessed May 16, 2025, https://www.researchgate.net/publication/352717865_A_Two-level_GREG_Estimator_for_Consistent_Estimation_in_Household_Surveys
24.	A Two‐level GREG Estimator for Consistent Estimation in Household Surveys - EconStor, accessed May 16, 2025, https://www.econstor.eu/bitstream/10419/284842/1/INSR_INSR12460.pdf
25.	Small Area Estimation: Informative Sampling and Two-fold Models - Carleton University Institutional Repository, accessed May 16, 2025, https://repository.library.carleton.ca/downloads/5712m731k?locale=en
26.	STATISTICAL MODELS FOR SMALL AREA ESTIMATION - Repositorio CIMAT, accessed May 16, 2025, https://cimat.repositorioinstitucional.mx/jspui/bitstream/1008/351/2/TE%20517.pdf
27.	A COMPREHENSIVE OVERVIEW OF UNIT-LEVEL MODELING OF SURVEY DATA FOR SMALL AREA ESTIMATION UNDER INFORMATIVE SAMPLING, accessed May 16, 2025, https://par.nsf.gov/servlets/purl/10447837
28.	Tree-Based Machine Learning in Small Area Estimation Patrick Krennmair 1, Nora W ¨urz 2 and Timo Schmid3, accessed May 16, 2025, https://isi-iass.org/home/wp-content/uploads/Survey_Statistician_2022_July_N86_03.pdf
29.	www150.statcan.gc.ca, accessed May 16, 2025, https://www150.statcan.gc.ca/n1/pub/12-001-x/2021002/article/00007-eng.pdf
30.	Fay–Herriot model - Wikipedia, accessed May 16, 2025, https://en.wikipedia.org/wiki/Fay%E2%80%93Herriot_model
31.	Small Area Estimation (Wiley Series in Survey Methodology): 9781118735787: Rao, J. N. K., Molina, Isabel - Amazon.com, accessed May 16, 2025, https://www.amazon.com/Small-Estimation-Wiley-Survey-Methodology/dp/1118735781
32.	Rao, J.N.K. and Molina, I. (2015) Small Area Estimation. 2nd Edition, Wiley, New York. - References - Scientific Research Publishing, accessed May 16, 2025, https://www.scirp.org/reference/referencespapers?referenceid=2066005
33.	Unifying area and unit-level small area estimation through calibration - arXiv, accessed May 16, 2025, https://arxiv.org/html/2403.15384v1
34.	(PDF) Model-Based Prediction of Finite Population Totals - ResearchGate, accessed May 16, 2025, https://www.researchgate.net/publication/251462184_Model-Based_Prediction_of_Finite_Population_Totals
35.	Generic small area estimation • SUMMER - GitHub Pages, accessed May 16, 2025, https://richardli.github.io/SUMMER/articles/web_only/small-area-estimation.html
36.	openknowledge.fao.org, accessed May 16, 2025, https://openknowledge.fao.org/server/api/core/bitstreams/b0a12b78-7263-45ab-931d-f27f772057b0/content
37.	Best linear unbiased prediction - Wikipedia, accessed May 16, 2025, https://en.wikipedia.org/wiki/Best_linear_unbiased_prediction
38.	Prediction with measurement errors in finite populations - PMC, accessed May 16, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC3230038/
39.	On Estimation and Prediction in Linear Mixed Models - DiVA portal, accessed May 16, 2025, https://www.diva-portal.org/smash/get/diva2:1826082/FULLTEXT01.pdf
40.	Small area benchmarked estimation under the basic unit level model when the sampling rates are non-negligible, accessed May 16, 2025, https://www150.statcan.gc.ca/n1/pub/12-001-x/2021001/article/00007-eng.pdf
41.	That BLUP is a Good Thing: The Estimation of Random Effects, accessed May 16, 2025, https://arnabc74.github.io/linmod/blup.pdf
42.	Calibration Estimation in Survey Sampling, accessed May 16, 2025, https://dr.lib.iastate.edu/bitstreams/3d17310a-c037-43e4-b683-b97493b8268b/download
43.	Calibration Estimators in Survey Sampling: Journal of the American ..., accessed May 16, 2025, https://www.tandfonline.com/doi/abs/10.1080/01621459.1992.10475217
44.	Calibration Estimators in Survey Sampling - Sci-Hub, accessed May 16, 2025, https://dacemirror.sci-hub.ru/journal-article/0de5e6874e4fdef5ab49807b2b862fe3/deville1992.pdf
45.	A Model-Calibration Approach to Using Complete Auxiliary ..., accessed May 16, 2025, https://www.researchgate.net/publication/4744252_A_Model-Calibration_Approach_to_Using_Complete_Auxiliary_Information_From_Survey_Data
46.	Neural Network Calibration Estimation with High Dimensional Data in Survey Sampling - Pakistan Journal of Statistics, accessed May 16, 2025, https://www.pakjs.com/wp-content/uploads/2023/07/39307.pdf
47.	Recent developments in calibration estimation - CiteSeerX, accessed May 16, 2025, https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=e901e5ababa0bebea8b19dd825b3f8be3cea0038
48.	Multivariate Regression Analysis | Stata Data Analysis Examples - OARC Stats - UCLA, accessed May 16, 2025, https://stats.oarc.ucla.edu/stata/dae/multivariate-regression-analysis/
49.	Gaussian Graphical Model Estimations in Multivariate Linear Regression: A method and applications in omics studies - European Society of Medicine, accessed May 16, 2025, https://esmed.org/MRA/mra/article/view/6088
50.	An alternative pseudolikelihood method for multivariate random-effects meta-analysis - PMC, accessed May 16, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC4305202/
51.	Empirical and Simulation-Based Comparison of Univariate and Multivariate Meta-Analysis for Binary Outcomes | Effective Health Care (EHC) Program, accessed May 16, 2025, https://effectivehealthcare.ahrq.gov/products/methods-binary-outcomes-evaluation/research
52.	MULTIVARIATE LOGISTIC-ASSISTED ESTIMATORS OF TOTALS FROM CLUSTERED SURVEY SAMPLES 1. INTRODUCTION - ResearchGate, accessed May 16, 2025, https://www.researchgate.net/profile/Richard-Valliant/publication/347982664_Multivariate_Logistic-Assisted_Estimators_of_Totals_From_Clustered_Survey_Samples/links/605de480299bf173676e7c4b/Multivariate-Logistic-Assisted-Estimators-of-Totals-From-Clustered-Survey-Samples.pdf
53.	S−estimators for functional principal component analysis - UBC Department of Statistics |, accessed May 16, 2025, https://www.stat.ubc.ca/~matias/Boente_Salibian_TR_06_06_2013.pdf
54.	Principal component analysis: a review and recent developments - PMC - PubMed Central, accessed May 16, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC4792409/