Weighted Likelihood Estimating Equations

• Weighted Likelihood Estimating Equations (WLEE) are score-based equations that replace the raw likelihood score with weighted contributions to enhance robustness.
• They adjust weights based on data-model compatibility, addressing issues like outliers, dependence structures, and class imbalances in a variety of applications.
• WLEE methods achieve asymptotic efficiency under correct model specification and serve as a versatile design pattern for robust and efficient inference.

Weighted Likelihood Estimating Equations (WLEE) are score-based estimating equations in which the ordinary likelihood score is replaced by a weighted score, typically to control robustness, encode design information, accommodate dependence, or optimize efficiency under a working model. In its most generic form, WLEE replaces the maximum-likelihood condition $U(\theta)=\sum_{i=1}^n s(X_i;\theta)=0$ by a weighted analogue in which each score contribution is multiplied by a scalar or matrix weight. Across the literature, this includes observation-wise robustification for i.i.d. models, matrix-weighted score equations for correlated data, power-likelihood reweighting in Bayesian prediction, calibration-based weights in empirical likelihood, and optimal estimating-function constructions for point processes (Majumder et al., 2016, Nikoloulopoulos, 2015, Agostinelli et al., 22 Jul 2025).

1. General formulation and scope

The canonical WLEE form in independent-data settings is

$$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$

where $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$ is the ordinary score and the weights $w_i$ may depend on the current parameter value, the observation, the empirical distribution, or external design information. In robust parametric work, the weights are data-adaptive and typically lie in $[0,1]$, with $w_i\approx 1$ when the observation is compatible with the model and smaller values when incompatibility is detected. In other strands, such as class-imbalance correction, the weights are exogenous and fixed in advance; in empirical-likelihood and survey problems they arise from constrained optimization; and in correlated-data models they are matrix-valued and derived from sensitivity and covariance calculations rather than from outlier diagnostics (Majumder et al., 2016, Lazic, 23 Apr 2025, Tan et al., 2020, Chaudhuri et al., 2022).

A second major formulation appears in dependent-data problems, where the weighting acts on stacked marginal scores rather than on scalar observation scores. For longitudinal GLM margins, the weighted scores method uses

$$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$

with $A_i^{(1)}$ the marginal sensitivity matrix and $\Sigma_{i,1}$ the working covariance of the univariate score vector under a proper multivariate working model. This construction is explicitly Godambe-oriented: the weights are chosen to precondition the score by its sensitivity and variability, rather than to suppress outliers (Nikoloulopoulos, 2015, Nikoloulopoulos, 2015).

These formulations are mathematically allied because each replaces the raw likelihood score by a weighted estimating function while retaining likelihood-based local geometry. What changes across applications is the source of the weights: discrepancy diagnostics, dependence structure, class proportions, calibration constraints, component densities, cellwise reliabilities, or predictable window functionals in a point process.

2. Robust score-weighting for independent data

A central line of development treats WLEE as a robust alternative to maximum likelihood. Majumder et al. propose

$$U_w(\theta)=\sum_{i=1}^n H(\tau_{n,\theta}(X_i))\,s(X_i;\theta)=0,$$

where the residual $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$0 compares the empirical distribution or survival function with the model only in the tails. Specifically, for a tail fraction $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$1, $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$2 equals $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$3 in the left tail, $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$4 in the right tail, and $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$5 in the central region, so central observations receive weight $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$6. The weight function $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$7 is required to satisfy $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$8, $$U_w(\theta)=\sum_{i=1}^n w_i(\theta;X_i)\,s(X_i;\theta)=0,$$9, $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$0, and $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$1, which yields full Fisher efficiency at the correctly specified model while selectively downweighting incompatible tail observations (Majumder et al., 2016).

The same paper emphasizes a notable subtlety: the first-order influence function of the resulting estimator coincides with that of the MLE, so robustness is not achieved through a bounded first-order IF. Instead, the robust behavior is explained by second-order bias reduction, together with the fact that only observations exhibiting empirical-model tail disagreement are downweighted. This distinguishes the method from schemes that continuously downweight a fixed fraction of observations regardless of compatibility (Majumder et al., 2016).

Other robust WLEE constructions use Pearson-type residuals based on smoothed densities or lower-dimensional summaries. Agostinelli and Greco replace full multivariate density comparison by comparison on the distribution of squared Mahalanobis distances, avoiding the curse of dimensionality that hampers density-based Pearson residuals in high dimension. Their multivariate location-scatter estimator solves weighted normal score equations in which the weights are built from a residual adjustment function applied to a univariate kernel estimate of the distance distribution under a $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$2 reference law; the resulting estimator is consistent and first-order efficient at the model, and it supports outlier detection and weighted PCA (Agostinelli et al., 2017).

A related depth-based line compares model depth and sample depth rather than densities. An earlier formulation defines

$s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$3

and plugs this quantity into a residual-to-weight map, with Tukey halfspace depth used in the multivariate normal examples (Agostinelli, 2018). A later formulation stabilizes the discrepancy through a Depth Pearson Residual,

$s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$4

for $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$5, studies asymptotic normality under the standard regularity conditions typically assumed for the MLE, proves full asymptotic efficiency at the model, and derives finite sample breakdown results for location and scatter in elliptically symmetric families. For the multivariate normal model, the half-space depth is explicitly linked to the squared Mahalanobis distance through $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$6 (Agostinelli et al., 22 Jul 2025).

3. Dependence-adjusted weighted scores for correlated data

In longitudinal discrete-data analysis, WLEE appears in a distinct but closely related form as a weighted scores method for GLM margins. For binary and count responses, the method starts from the stacked marginal score vector $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$7 and weights it using a working discretized multivariate normal (DMVN) model. The DMVN supplies a proper multivariate model for discrete outcomes through multivariate normal rectangle probabilities, with exchangeable, AR(1), and unstructured latent correlations available. The resulting estimating equation,

$s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$8

is interpreted as a Godambe-optimal preconditioning of the independent score among estimating functions built from the marginal scores (Nikoloulopoulos, 2015).

This framework is positioned as a competitor to generalized estimating equations. Both methods remain consistent for $s(X_i;\theta)=\nabla_\theta \log f(X_i;\theta)$9 when the marginal GLM is correctly specified and clusters are independent, but the weighted scores method uses a proper multivariate working model, so its dependence parameters are interpretable latent correlations and remain within feasible bounds. By contrast, a GEE working correlation for discrete data need not correspond to any proper joint model and can violate Fréchet bounds. The paper further couples WLEE with a pairwise composite-likelihood step, denoted CL1, to estimate the working dependence and to define composite-likelihood information criteria

$w_i$0

$w_i$1

which are used for both correlation-structure selection and variable selection. In the reported simulations, CL1-based criteria outperform GEE-based criteria overall; in Hamilton’s depression data, they favor an AR(1) structure and show that misspecifying the working structure in GEE can alter inferential conclusions (Nikoloulopoulos, 2015).

The same weighted-scores philosophy was extended to longitudinal ordinal data. There the key advantage is structural: the method works directly with ordinal margins and their univariate scores, rather than first expanding each ordinal response into $w_i$2 binary indicators as in standard ordinal GEE. This avoids the large $w_i$3-dimensional covariance matrices and the need to specify cross-category indicator correlations. The working weight remains of the form $w_i$4, now built from ordinal univariate information and DMVN-based bivariate rectangle probabilities. The reported timing studies show substantial computational gains over indicator-based GEE, with weighted scores 9–72× faster in large-$w_i$5 or large-$w_i$6 settings, while retaining near-ML efficiency and robust type-I error under dependence misspecification (Nikoloulopoulos, 2015).

4. Calibration, design, and purpose-specific weighting

Not all WLEE are motivated by robustness. In Bayesian prediction under class imbalance, weighted likelihood is introduced as a power-likelihood in which each observation contributes $w_i$7, so that

$w_i$8

For inverse-proportion weighting, if $w_i$9 is the empirical class proportion of class $[0,1]$0, a convenient normalized choice is $[0,1]$1. In logistic, multinomial softmax, and ordered logistic models, this directly rescales the score and Hessian. The reported experiments show unchanged AUC in one binary example but improved balanced accuracy, sensitivity, F1, and P4, together with worse Brier score and worse calibration, making the trade-off explicitly application-dependent (Lazic, 23 Apr 2025).

Empirical-likelihood approaches generate WLEE through constrained optimization rather than by specifying a weight function directly. For average treatment effects in randomized trials, two-sample empirical likelihood imposes covariate-balance constraints separately in treated and control groups, yielding weights of the form

$[0,1]$2

and the weighted ATE estimator is the difference of the two weighted sample means. When the working outcome regressions are correctly specified, the resulting estimator attains the semiparametric efficiency bound; under missing at random it can be doubly robust, and with multiple working models it can be multiply robust (Tan et al., 2020).

A related conditional empirical-likelihood construction appears for informative complex surveys. There, the conditional inclusion probability $[0,1]$3 leads to WLEE of the form

$[0,1]$4

when no auxiliary calibration constraint is imposed, and to

$[0,1]$5

when population-level estimating equations are added. This provides a likelihood-based derivation of inverse–conditional-probability weighting and its calibration-adjusted version, and the paper shows asymptotic efficiency gains relative to other probability-weighted analogs (Chaudhuri et al., 2022).

Other purpose-specific reweightings extend the same logic. In cellwise weighted data, the likelihood is defined by “unpacking” each row into partially observed rows with row weights $[0,1]$6, reducing cellwise weighting to standard rowwise weighted observed likelihood on an enlarged incomplete-data set. For multivariate normal models this yields an EM algorithm and fast explicit approximations, denoted cwMean and cwCov, that are asymptotically equivalent to the cellwise maximum likelihood estimator (Rousseeuw, 2022). In robust mixture fitting, weighted complete estimating equations use component-specific weights $[0,1]$7 together with bias-correction terms $[0,1]$8 and $[0,1]$9, leading to an expectation–estimating-equation algorithm that robustifies Gaussian mixtures, mixtures of experts, and skew-normal mixtures (Sugasawa et al., 2020).

5. Asymptotic theory, optimal weighting, and computation

Across these formulations, the dominant asymptotic language is the Godambe or sandwich calculus. For an estimating function $w_i\approx 1$0, one defines

$w_i\approx 1$1

Then $w_i\approx 1$2, with plug-in sandwich estimators replacing $w_i\approx 1$3 and $w_i\approx 1$4 by empirical counterparts. This appears explicitly in longitudinal weighted scores, where the asymptotic covariance of $w_i\approx 1$5 is obtained from the Godambe information of the weighted score, and in the more classical robust WLEE literature, where asymptotic normality follows from Taylor expansion of the weighted score around the truth (Nikoloulopoulos, 2015, Majumder et al., 2016).

At the correctly specified model, several constructions recover full MLE efficiency by design. Majumder et al. obtain $w_i\approx 1$6 and $w_i\approx 1$7 because $w_i\approx 1$8 and $w_i\approx 1$9, so $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$0 converges to $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$1. The stabilized depth-based estimator proves the same limit under standard regularity assumptions and mild depth-process conditions, including settings where the data dimension diverges with the sample size. In both cases the weights are asymptotically inactive at the model and only alter behavior under contamination or misspecification (Majumder et al., 2016, Agostinelli et al., 22 Jul 2025).

The most explicit optimal-weighting result appears for compact-memory Hawkes processes. There, every sufficiently regular predictable functional of the recent lag window generates an unbiased compensator-based estimating equation

$$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$2

The likelihood score itself is one such member, with weight $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$3, and it provides the efficient benchmark. For a finite library $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$4, the paper defines

$$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$5

so the optimal asymptotic covariance is $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$6. A feasible two-step procedure estimates $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$7 from a first-step fit and then re-solves the estimating equation with $$\sum_{i=1}^n X_i^\top W_{i,\mathrm{work}}^{-1}s_i^{(1)}(\beta)=0, \qquad W_{i,\mathrm{work}}^{-1}=A_i^{(1)}(\beta)\Sigma_{i,1}(\beta,R)^{-1},$$8, exactly paralleling optimal GMM and optimal GEE weighting. The paper also gives a projection identity that expresses the efficiency loss of a library as the component of the score outside the predictable span generated by that library (Davis et al., 22 Jun 2026).

Computation is correspondingly heterogeneous. Robust parametric WLEE typically use Newton–Raphson or Fisher-scoring updates together with bootstrap root search when multiple solutions are possible. Longitudinal weighted scores require repeated evaluation of bivariate MVN rectangle probabilities and, for some information criteria, trivariate or four-variate margins. Cellwise weighting relies on EM after unpacking, mixture-model WCE on EEE iterations, and depth-based WLEE on repeated evaluation of halfspace depth or related depth functionals. The computational burden is therefore usually concentrated not in the score itself, but in the construction of the weights.

6. Applications, limitations, and interpretive issues

WLEE have been deployed in a wide range of problem classes. Longitudinal binary, count, and ordinal models use them to replace or extend GEE while retaining population-averaged interpretation and enabling likelihood-grounded selection of working correlation and covariates (Nikoloulopoulos, 2015, Nikoloulopoulos, 2015). Insurance-loss GLMs use score-based weighted likelihood to obtain robust estimation under heavy-tailed contamination, and extend the same mechanism to censored and truncated losses through weighted scores involving truncated densities and censoring-interval probabilities (Fung, 2022). Wrapped models for multivariate circular data use torus-based or unwrapped Pearson-type residuals, showing that WLEE can be built either directly on the torus or on an unwrapped Euclidean representation; the paper reports similar robust fits for both, with unwrapped implementations substantially faster (Agostinelli et al., 2024).

Several misconceptions recur. One is that WLEE necessarily means outlier downweighting. The literature considered here shows a broader picture: some WLEE are robust, but others are efficiency-oriented matrix weights for dependent data, fixed class weights for predictive imbalance, calibration weights for causal or survey inference, or component-wise weights for latent-structure models. Another misconception is that robustness in WLEE always corresponds to a bounded first-order influence function. Majumder et al. explicitly show that their estimator shares the MLE’s first-order IF at the model, with robustness emerging through second-order behavior and selective tail downweighting rather than first-order boundedness (Majumder et al., 2016).

The principal limitations are also recurrent. Many robust WLEE can have multiple roots, making initialization and root selection decisive; this is documented for tail-ratio weighting, depth-based weighting, and mixture-model WCE. Methods based on smoothed residuals or depths inherit tuning sensitivity through bandwidths, cutoffs, or smoothing kernels. Dependence-adjusted weighted scores rely on a working multivariate model; if the marginal mean model is misspecified, consistency is lost just as in other semiparametric estimating-equation methods. Computational cost can be substantial when the weights require MVN rectangle probabilities, multivariate depth computation, or repeated kernel estimates. Some explicit approximations, such as cwCov in cellwise weighting, need not be positive semidefinite in finite samples and may require regularization (Agostinelli et al., 2017, Rousseeuw, 2022, Sugasawa et al., 2020).

A plausible general implication is that WLEE are best viewed not as a single estimator but as a design pattern: retain the score structure of likelihood inference, then alter the effective contribution of data points, score components, or clusters by weights chosen for robustness, calibration, dependence adjustment, or efficiency. What unifies the family is not a specific weight formula, but the replacement of the raw score by a weighted estimating function whose asymptotic behavior is typically analyzed through Godambe information and whose practical success depends critically on how the weights encode model-data compatibility or inferential objectives.