Fact-checked by Grok 2 weeks ago

Weighted least squares

Weighted least squares (WLS) is a of ordinary least squares (OLS) that assigns different weights to observations based on their estimated variances, allowing for the analysis of with heteroscedastic errors where the variance of residuals is not constant across observations. This method minimizes the weighted sum of squared residuals, \sum w_i (y_i - \hat{y}_i)^2, where weights w_i are typically set to the inverse of the error variance for each point, w_i = 1/\sigma_i^2, thereby giving greater influence to more precise observations and reducing the impact of those with higher variability. In contrast to OLS, which assumes homoscedasticity and treats all observations equally, WLS addresses violations of this assumption by incorporating a diagonal weight matrix W into the estimation process, yielding the parameter \hat{\beta}_{WLS} = (X^T W X)^{-1} X^T W y. Under the assumptions of the Gauss-Markov theorem—linearity, unbiasedness, and known error covariance structure—the WLS is the best linear unbiased (BLUE), achieving minimum variance among all linear unbiased estimators. WLS is particularly valuable in applications such as , where error variances may vary with the scale of predictors, and in for analyzing grouped or clustered data with differing sample sizes. When the true variances are unknown, they can be estimated from residuals or auxiliary models, leading to feasible as an extension for correlated errors. This flexibility makes WLS a robust tool for improving the and reliability of estimates in real-world datasets exhibiting non-constant error structures.

Fundamentals

Definition and Relation to Ordinary Least Squares

Weighted least squares (WLS) is a generalization of (OLS) estimation in , designed to minimize a weighted sum of squared residuals when observations have unequal variances or differing levels of importance. In this approach, each data point is assigned a positive weight w_i, typically proportional to the inverse of its error variance, to give greater to more reliable observations. The core distinction from OLS lies in the objective function: OLS minimizes the unweighted sum \sum_{i=1}^n (y_i - \hat{y}_i)^2, assuming homoscedastic errors with equal weight 1 for all residuals, whereas WLS minimizes \sum_{i=1}^n w_i (y_i - \hat{y}_i)^2 to address heteroscedasticity or unequal precision. This weighting ensures that the estimator accounts for varying reliability across data points, leading to more efficient parameter estimates under non-constant variance conditions. WLS operates within the framework of the y = X\beta + \epsilon, where y is the n \times 1 response vector, X is the n \times p , \beta is the p \times 1 parameter vector, and \epsilon is the error vector with mean zero and diagonal \Sigma = \diag(\sigma_1^2, \dots, \sigma_n^2), reflecting independent but heteroscedastic errors. The method has roots in 19th-century astronomy, where introduced weighted least squares in 1821–1823 to refine estimates from observations of varying precision in and . It was later formalized in statistical contexts during the early , notably by William G. Cochran in 1937, who addressed weighting for combining results from similar experiments with heterogeneous variances.

Assumptions and Prerequisites

To apply weighted least squares (WLS) effectively, readers should possess foundational knowledge of linear regression models, including matrix notation such as the design matrix \mathbf{X} (an n \times p matrix of predictors) and the response vector \mathbf{y} (an n \times 1 vector of observations). This prerequisite enables understanding of how WLS extends ordinary least squares (OLS) by incorporating a weight matrix into the estimation process. The core assumptions of WLS mirror those of OLS in several respects but adapt the error variance structure to account for heteroscedasticity. Specifically, the model assumes linearity in the parameters, meaning the E(\mathbf{y} | \mathbf{X}) = \mathbf{X} \boldsymbol{\beta}, where \boldsymbol{\beta} is the . It also requires no perfect among the columns of \mathbf{X}, ensuring the weighted is invertible and estimates are uniquely defined. The errors are assumed to be independent across observations, with E(\boldsymbol{\epsilon}) = \mathbf{0}, though extensions like can accommodate specified correlation structures. Crucially, the variance of the errors is heteroscedastic, modeled as \text{Var}(\mathbf{y} | \mathbf{X}) = \boldsymbol{\Sigma}, a where the i-th diagonal element is \sigma_i^2; the weights w_i are then proportional to the inverse variances, w_i = 1 / \sigma_i^2, reflecting the relative of each observation. These weights must be known or reliably estimated in advance. In contrast to OLS, which assumes homoscedasticity (\sigma_i^2 = \sigma^2 for all i), WLS relaxes this by allowing varying error variances while retaining unbiasedness under the correct linear model. However, if homoscedasticity actually holds, applying WLS with non-uniform weights leads to inefficient estimates compared to OLS, though the point estimates remain unbiased. Misspecifying the weights—such as using incorrect inverse variances—does not introduce bias in the parameter estimates, provided the linearity assumption is met, but it results in a loss of efficiency and invalid standard errors, compromising inference. Violations of linearity or perfect multicollinearity can produce biased or unstable estimates, underscoring the need to verify these conditions prior to WLS application.

Motivation

Handling Heteroscedasticity

Heteroscedasticity refers to the situation in models where the variance of the error terms, denoted as \sigma_i^2, is not constant across observations but instead varies systematically, often depending on the values of the predictors or other factors. This violates a key assumption of (OLS), which requires homoscedasticity for optimal properties. Under heteroscedasticity, observations with larger error variances contribute disproportionately to the estimation, leading to inefficient OLS estimates despite their remaining unbiased. The presence of heteroscedasticity has significant consequences for in OLS . While the point estimates of the are unbiased, their variances are inflated, making OLS inefficient compared to alternatives that account for varying variances. More critically, the conventional errors computed under OLS are biased, often underestimated, which invalidates t-tests, F-tests, and intervals, potentially leading to incorrect conclusions about significance. Weighted least squares (WLS) mitigates this by assigning weights inversely proportional to the variances, transforming the model to restore homoscedasticity and yielding the best linear unbiased estimator () when the weights are correctly specified, as established by Aitken's of the Gauss-Markov theorem. Detecting heteroscedasticity is essential before applying WLS. One common approach is the Breusch-Pagan test, which involves regressing the squared OLS residuals on the independent variables and testing whether the coefficients are jointly zero using a ; a significant result indicates variance dependence on predictors. Another is White's test, which extends this by including squared and cross-product terms of the predictors in the auxiliary regression to capture more general, unspecified forms of heteroscedasticity, also evaluated via a . These tests help confirm the violation without assuming a specific variance structure. A classic real-world illustration of heteroscedasticity appears in regressions modeling household expenditure as a function of . For lower- households, food spending tends to closely around predicted values due to constraints, resulting in low variance. In contrast, higher- households exhibit greater variability in expenditures, influenced by diverse choices in dining and grocery habits, leading to increasing \sigma_i^2 with income levels. Applying WLS with weights based on this pattern, such as inversely related to income, can then produce more efficient estimates and reliable inference for such data.

Incorporating Prior Knowledge via Weights

In weighted least squares (WLS), weights can extend beyond correcting for heteroscedasticity to incorporate non-variance-based prior knowledge, such as the relative reliability of observations, effective sample sizes in aggregated , or subjective assessments from experts. For instance, in analyses of where individual are aggregated into categories (e.g., brackets or groups), weights are often set proportional to the number of underlying in each group, reflecting the greater informational or of larger aggregates. This approach ensures that estimates prioritize groups with more robust empirical support, drawing on the understanding that variance decreases with sample size under standard assumptions. Similarly, weights may be chosen inversely proportional to the of obtaining an or the estimated in auxiliary sources, allowing practitioners to emphasize more economical or accurate inputs in resource-constrained settings. A prominent application of such non-variance weighting occurs in , where WLS combines effect estimates from multiple studies by assigning weights based on study-specific sample sizes rather than solely on inverse variances. When individual study variances are similar or difficult to estimate precisely, weighting by sample size provides a practical that favors larger, more precise studies, leading to an overall estimate that better reflects the accumulated evidence across the literature. This method has been shown to involve a compared to , being less biased but less efficient, particularly when sample sizes serve as reliable indicators of precision. The primary advantage of incorporating prior knowledge through these weights lies in its ability to enhance estimator efficiency by integrating domain-specific insights, even in scenarios of homoscedasticity where ordinary least squares would suffice. By leveraging expert judgment or structural information—such as incomplete priors on parameter relations—WLS can yield more informed predictions than unweighted methods, as demonstrated in mixed estimation frameworks that blend sample data with subjective beliefs. However, a key limitation is that subjective or auxiliary-based weights, if misspecified, can distort the estimator's properties; while unbiasedness is generally preserved under correct model specification, suboptimal weights may inflate variance, reducing efficiency. This underscores the need for careful validation of weight choices to avoid compromising inferential reliability.

Mathematical Formulation

Model Specification

The weighted least squares (WLS) model is formulated within the framework of as y = X\beta + [\epsilon](/page/Epsilon), where y is an n \times 1 of observed responses, X is an n \times p full-rank with rows corresponding to covariates for each observation, \beta is a p \times 1 of unknown coefficients, and [\epsilon](/page/Epsilon) is an n \times 1 error satisfying E([\epsilon](/page/Epsilon)) = 0 and \text{Var}([\epsilon](/page/Epsilon)) = \Sigma, with \Sigma being a positive definite n \times n . This setup generalizes ordinary least squares by allowing for non-constant error variances or correlations, where \Sigma captures heteroscedasticity or dependence in the errors. The core of WLS estimation involves minimizing the weighted sum of squared residuals, expressed as the objective function (y - X\beta)^T W (y - X\beta), where W = \Sigma^{-1} is the n \times n weight matrix that inversely scales the residuals according to their error variances; when errors are uncorrelated, W is diagonal with entries w_i = 1 / \sigma_i^2 for the i-th observation's variance \sigma_i^2. This weighting ensures that observations with lower variance contribute more to the parameter estimates, aligning with the Gauss-Markov theorem for efficiency under known \Sigma. Equivalently, the WLS problem can be transformed into an ordinary least squares (OLS) estimation on adjusted variables by defining y^* = W^{1/2} y and X^* = W^{1/2} X, where W^{1/2} is the symmetric positive definite square root of W; the objective then reduces to minimizing \| y^* - X^* \beta \|^2, which preserves the linear structure while accounting for the error covariance. This transformation highlights how WLS reweights the data to achieve homoscedasticity in the transformed space.

Weight Selection and Matrix Construction

In weighted least squares (WLS), the ideal weights for observations with uncorrelated errors are chosen as the reciprocals of the error variances, w_i = 1/\sigma_i^2, where \sigma_i^2 represents the variance of the error term for the i-th observation; this assignment ensures that more precise observations (lower variance) receive higher influence in the estimation process. The weight matrix W is then constructed as a diagonal matrix with these w_i on the diagonal, W = \operatorname{diag}(w_1, w_2, \dots, w_n), which simplifies computation while assuming independence among errors. When errors are correlated, the weight matrix extends to the inverse of the full of the errors, W = \Sigma^{-1}, transforming WLS into (GLS); this full matrix accounts for both heteroscedasticity and serial or spatial correlations, though its estimation is more complex and typically deferred to GLS-specific methods. For practical implementation, the diagonal form of W is preferred in standard WLS due to its computational efficiency and sufficiency under the uncorrelated errors assumption. Since error variances are often unknown in practice, weights must be estimated iteratively or through modeling. A common approach begins with an initial ordinary least squares (OLS) fit to compute residuals e_i = y_i - \hat{y}_i, from which preliminary variance estimates \hat{\sigma}_i^2 are derived, such as by squaring the residuals or regressing them against predictors or fitted values to model the variance function. These estimates yield initial weights \hat{w}_i = 1/\hat{\sigma}_i^2, which are then used in a WLS fit; the process iterates by updating residuals and weights until convergence, often requiring only one or two iterations for stability. Alternatively, model-based estimation assumes a form for the variance, such as \sigma_i^2 = \sigma^2 h(x_i) where h(\cdot) is a known function (e.g., h(x_i) = 1 + 0.5 x_i^2), allowing weights to be constructed directly from the predictors. Common pitfalls in weight selection include the variance model from limited data, which can lead to unstable or biased estimates, and failing to transform variables appropriately, resulting in an invalid positive-definite W. Small sample sizes exacerbate these issues, as estimated weights become sensitive to outliers or minor variations in residuals, potentially undermining the efficiency gains of WLS.

Estimation Methods

Closed-Form Solution

The closed-form solution for the weighted least squares (WLS) estimator \hat{\beta} is derived by minimizing the objective function, which is the weighted sum of squared residuals S(\beta) = (y - X\beta)^T W (y - X\beta), where y is the response vector, X is the , \beta is the parameter vector, and W is the positive definite diagonal weight matrix. To obtain the minimizer, differentiate S(\beta) with respect to \beta and set the result to zero: \frac{\partial S}{\partial \beta} = -2 X^T W (y - X \beta) = 0. This yields the normal equations: X^T W X \beta = X^T W y. Solving for \beta gives the explicit estimator: \hat{\beta} = (X^T W X)^{-1} X^T W y. This solution exists provided that X^T W X is invertible, which requires the design matrix X to have full column rank and the weight matrix W to be positive definite. For numerical computation, direct matrix inversion of X^T W X can be unstable due to potential ill-conditioning; instead, a QR decomposition of the transformed W^{1/2} X is recommended, as it provides a stable alternative by reducing the problem to an solution on the transformed .

Iterative and Numerical Approaches

In cases where the weight W depends on the unknown \beta, such as in variants, the closed-form solution for weighted least squares (WLS) becomes infeasible, necessitating iterative methods to achieve convergence. (IRLS) addresses this by alternating between estimating \beta using WLS with current weights and updating the weights based on the residuals from that estimate, repeating until the parameter values stabilize within a specified tolerance. This approach, originally developed for to downweight outliers via functions like Huber's, ensures that weights reflect the reliability of each observation iteratively. Numerical challenges arise in WLS computations, particularly when the matrix X^T W X is ill-conditioned due to in X or extreme variation in the diagonal elements of W, leading to unstable or imprecise parameter estimates. To mitigate this, regularization techniques such as can be incorporated by adding a penalty term to the objective function, or the Moore-Penrose pseudoinverse can be employed to compute a least-squares solution that handles rank deficiency. These methods stabilize the inversion process without assuming full rank, though they may introduce that requires careful tuning. Software packages facilitate efficient WLS implementation, including iterative variants like IRLS. In , the lm function supports weights via the weights argument, performing WLS directly or iteratively if needed through extensions like MASS::rlm for robust cases. Similarly, Python's statsmodels provides the WLS class in statsmodels.regression.linear_model, which accepts a weights and handles internally via optimized linear algebra routines. IRLS and related numerical approaches are particularly suited for large datasets where direct matrix inversion is computationally prohibitive, scenarios with non-diagonal covariance structures approximated by diagonal weights, or when weights must be refined based on provisional \beta estimates, as in robust or heteroscedasticity-adjusted models.

Statistical Properties

Estimator Variance and Covariance

The covariance matrix of the weighted least squares (WLS) estimator \hat{\beta} is given by \operatorname{Var}(\hat{\beta}) = \sigma^2 (X^T W X)^{-1}, where \sigma^2 is the error variance, X is the , and W is the diagonal weight matrix, assuming the weights are correctly specified to account for the error structure. This formula arises from the linear form of the estimator \hat{\beta} = (X^T W X)^{-1} X^T W y and the assumption that the errors have mean zero and \sigma^2 W^{-1}. The diagonal elements of \operatorname{Var}(\hat{\beta}) represent the variances of the individual parameter estimates \hat{\beta}_j, quantifying the uncertainty in each coefficient, while the off-diagonal elements capture the covariances between pairs of estimates, indicating how the estimation of one parameter affects others due to shared data structure. These covariances are particularly useful for understanding multicollinearity or dependencies in the parameter space. When the weight matrix W = \Sigma^{-1}, where \Sigma is the true covariance matrix, the WLS achieves the minimum variance among all linear unbiased s, as established by the Gauss-Markov , and in the Gaussian case, it attains the Cramér-Rao lower bound for efficiency. Under heteroscedasticity, this results in smaller variances compared to least squares (OLS) estimates, which ignore varying variances and thus fail to minimize risk in such scenarios. Additionally, under standard regularity conditions such as fixed , bounded moments, and correct weight specification, the WLS is asymptotically normal: \sqrt{n} (\hat{\beta} - \beta) \xrightarrow{d} N(0, \lim_{n \to \infty} n \sigma^2 (X^T W X / n)^{-1}). The error variance \sigma^2 is typically estimated from the weighted residuals as \hat{\sigma}^2 = \frac{(y - X \hat{\beta})^T W (y - X \hat{\beta})}{n - p}, where n is the sample size and p is the number of parameters; this unbiased estimator adjusts the weighted sum of squared residuals by the degrees of freedom to provide a consistent measure of dispersion. Substituting \hat{\sigma}^2 into the covariance matrix yields the estimated \widehat{\operatorname{Var}}(\hat{\beta}), which is used for inference on the parameters.

Confidence Intervals for Parameters

In weighted least squares (WLS), inference for the parameters \hat{\beta} relies on the variance-covariance matrix of the estimator, which provides the foundation for standard errors and subsequent confidence intervals. The standard error for the j-th parameter estimate is given by \text{SE}(\hat{\beta}_j) = \sqrt{ \left[ (X^T W X)^{-1} \right]_{jj} \hat{\sigma}^2 }, where \hat{\sigma}^2 is the estimated error variance, typically computed as the weighted mean squared error from the residuals. This measure quantifies the precision of \hat{\beta}_j after accounting for the heteroscedasticity addressed by the weight matrix W. For finite samples, confidence intervals for individual parameters assume approximate of the and use the t-distribution: \hat{\beta}_j \pm t_{n-p, 1-\alpha/2} \cdot \text{SE}(\hat{\beta}_j), with n - p , where n is the number of observations and p is the number of parameters. This construction ensures that the interval captures the true \beta_j with probability $1 - \alpha, adjusting for the reduced in the WLS framework. Hypothesis testing for parameters in WLS follows standard procedures adapted to the weighted model. For a single coefficient, the Wald test statistic \hat{\beta}_j / \text{SE}(\hat{\beta}_j) follows a t-distribution with n - p degrees of freedom under the null hypothesis H_0: \beta_j = 0. For joint tests on subsets of \beta involving q restrictions, an F-test is used, with the statistic distributed as F_{q, n-p} under the null, comparing restricted and unrestricted weighted residual sums of squares. When weights are estimated rather than known—such as in feasible —inference based on the conventional standard errors may be invalid due to uncertainty in W. In these cases, robust standard errors via the sandwich estimator provide consistent inference, adjusting the to account for potential misspecification: \widehat{\text{Var}}(\hat{\beta}) = (X^T W X)^{-1} (X^T W \hat{\Omega} W X) (X^T W X)^{-1}, where \hat{\Omega} is a of squared weighted residuals (often with finite-sample like HC3). This approach maintains valid confidence intervals and tests even under weight estimation error.

Residual Diagnostics

In weighted least squares (WLS) regression, residual diagnostics begin with the computation of weighted residuals, defined as r_i = \sqrt{w_i} (y_i - \hat{y}_i), where w_i is the weight for the i-th observation, y_i is the observed response, and \hat{y}_i is the fitted value. These weighted residuals standardize the ordinary residuals by the square root of the weights, accounting for the heteroscedasticity assumed in the model. If the weighting scheme correctly captures the variance structure, the weighted residuals should exhibit homoscedasticity (constant variance) and approximate a normal distribution with mean zero. Common diagnostic tools include graphical assessments and formal statistical tests applied to these weighted residuals. Plots of weighted residuals against fitted values or predictor variables should display no discernible patterns, such as funnel shapes indicative of unresolved heteroscedasticity, to confirm model adequacy. For , a Q-Q plot of the weighted residuals versus theoretical quantiles can visually check for deviations from a straight line, while the Shapiro-Wilk test provides a formal test; under the of , the W approaches 1, with small p-values signaling non-normality. These diagnostics help verify that the assumptions underlying WLS hold, ensuring reliable . Influence measures in WLS adapt standard regression diagnostics to incorporate weights, identifying observations that disproportionately affect parameter estimates. The weighted version of Cook's distance, D_i = \frac{ r_i^{*2} h_{ii} }{ p (1 - h_{ii}) }, where r_i^* is the studentized weighted residual and h_{ii} is the i-th leverage adjusted for weights, quantifies the change in fitted values across all observations when the i-th case is excluded; values exceeding $4/n (with n observations and p parameters) flag influential points. Similarly, weighted DFFITS, which measures the scaled difference in predicted values with and without the observation, uses the formula \text{DFFITS}_i = r_i^* \sqrt{h_{ii}/(1 - h_{ii})}, with thresholds around $2\sqrt{p/n} indicating influence. These metrics, computed via weighted leverage matrices, prioritize observations with higher weights in assessing impact. To validate the model, examine whether the weighting has successfully addressed heteroscedasticity by re-plotting weighted residuals against fitted values or predictors; random scatter around zero confirms resolution, while persistent trends suggest misspecified weights. If diagnostics reveal remaining issues, such as patterns implying correlated errors beyond simple heteroscedasticity, generalized least squares (GLS) may be appropriate as an extension, incorporating a full structure for the errors.

Applications and Extensions

Practical Examples

One practical application of weighted least squares (WLS) arises in linear regression scenarios with heteroscedastic errors, where the variance of the residuals varies across observations. A classic example is Francis Galton's 1877 study on the inheritance of pea plant characteristics, using data on parent and progeny pea diameters to model the relationship between parental size and average offspring size. In this dataset of seven observations, the standard deviation (SD) of the progeny measurements increases with parental size, indicating heteroscedasticity. To apply WLS, weights are computed as the inverse of the squared standard deviations (w_i = 1 / \sigma_i^2), where \sigma_i is the for each parental group; this downweights observations with higher variability and upweights those with lower variability. The WLS model is then estimated by minimizing the weighted sum of squared residuals, yielding the fitted \hat{Y} = 0.12796 + 0.2048 X, where Y is progeny and X is parent . In comparison, ordinary least squares (OLS) produces \hat{Y} = 0.12703 + 0.2100 X, which overemphasizes noisier right-side points. The WLS fit pulls the line closer to low-variance points on the left, reducing the overall variance and improving without biasing the coefficients. This adjustment enhances predictive accuracy; for instance, WLS better captures the inheritance pattern by accounting for measurement precision, leading to more reliable inferences about genetic . Another common use of WLS is in analyzing survey data to adjust for complex sampling designs, such as certain subgroups, through post-stratification weights that align sample estimates with known totals. Consider the 2009-2010 and Nutrition Examination Survey (NHANES), a large-scale U.S. health study with 10,162 respondents, where sampling weights w_i correct for unequal selection probabilities (e.g., of children and minorities). These weights are normalized so their sum equals the sample size, and for WLS, the diagonal weight matrix uses the positive square roots \sqrt{w_i} to incorporate design effects into the estimation. The model regresses body weight (kg) on (cm), with WLS solving the weighted normal equations (X^T W X) \hat{\beta} = X^T W y, where W is diagonal with \sqrt{w_i}. Using \sqrt{w_i}, WLS yields a population-representative body weight of 70.86 kg and, in a no-intercept model, a of 0.46, whereas unweighted gives a biased mean of 63.29 kg (underestimating due to oversampling lighter subgroups) and a slope of 0.47. WLS thus produces unbiased, design-consistent estimates, changing coefficients to better represent the broader population and improving predictions for public health applications like obesity modeling. In both cases, WLS refines coefficient estimates and variance properties, yielding more precise predictions tailored to and sampling structure compared to unweighted methods.

Relation to Generalized Least Squares

(GLS) provides a broader framework for that accommodates both heteroscedasticity and in the error terms, positioning weighted least squares (WLS) as a specific instance within it. The GLS minimizes the weighted sum of squared residuals given by (y - X\beta)^T \Sigma^{-1} (y - X\beta), where \Sigma denotes the full n \times n covariance matrix of the error vector, assuming \mathbb{E}(\epsilon) = 0 and \text{Var}(\epsilon) = \Sigma. This approach yields the BLUE (best linear unbiased estimator) under the Gauss-Markov theorem when \Sigma is known, enhancing efficiency over ordinary least squares by accounting for the error structure. WLS emerges as a special case of GLS when \Sigma is diagonal, implying uncorrelated errors with potentially varying variances (heteroscedasticity) but no off-diagonal covariances. In this scenario, the inverse weights are the diagonal elements of \Sigma^{-1}, reducing GLS to the familiar WLS form where each observation is scaled by its variance. GLS is particularly advantageous over WLS in settings with correlated errors, such as in data or intra-cluster correlations in panel or grouped observations, where WLS would fail to capture dependencies and lead to inefficient estimates. For instance, in autoregressive models, GLS transforms the data to whiten the errors, yielding more precise estimates compared to WLS, which assumes . When \Sigma is unknown, feasible GLS (FGLS) addresses this by iteratively estimating the covariance structure—often starting with an initial ordinary least squares fit to obtain residuals, then updating weights via maximum likelihood or method of moments—and re-estimating the model until , akin to the procedure. This approximation performs nearly as well as true GLS under mild conditions and is widely implemented in statistical software for practical analysis. Historically, WLS preceded the GLS formulation; A. C. Aitken introduced the weighted least squares estimator in 1936 as the minimizer of a weighted sum of squared residuals for linear combinations of observations with known variances. The generalization to arbitrary covariance structures was advanced by in 1965, who developed a unified theory for models with parameters, encompassing GLS as a cornerstone for inference in correlated error settings.