Leverage Contamination in Robust Regression Inference

Why Leverage Changes the Meaning of Robustness

A robust regression procedure may limit the effect of a response outlier and still be vulnerable to a case whose predictors are extreme. These are different forms of contamination. A response outlier has an unusual response relative to its predictors; a high-leverage case has an unusual predictor vector and can exert substantial influence even when the residual score is bounded. The distinction matters for inference: a robust coefficient estimate and a sandwich standard error do not, by themselves, guarantee a confidence region for the uncontaminated coefficient.

Write \(Y_i\in\mathbb{R}^p\) for a response vector and \(x_i\in\mathbb{R}^q\) for a predictor vector, including an intercept if required. The regression coefficient matrix is \(B\in\mathbb{R}^{q\times p}\), with residual \(e_i(B)=Y_i-B^Tx_i\). The target in a model-based analysis is often a specified coefficient \(B_0\); under contamination or misspecification, an estimating procedure may instead converge to a different population value.

Definition: A response outlier is an observation with a large residual relative to the fitted regression relation. A high-leverage case has a predictor vector unusual relative to the predictor distribution or design. Leverage is a property of the predictors, not a synonym for a large residual.

Estimating Scores and Leverage

A common robust regression equation has the form

\sum_{i=1}^n g_i(B)=0,\qquad g_i(B)=w(x_i)\{\psi(e_i(B))\otimes x_i\}.

Here \(\psi:\mathbb{R}^p\to\mathbb{R}^p\) is a residual score, \(w(x_i)\) is an optional predictor weight, and \(\otimes\) denotes the Kronecker product. A radial Huber score, for example, is \(\psi(e)=e\min(1,c/\|e\|)\), with \(\psi(0)=0\). It is bounded in norm by \(c\). In the scalar-response case, the score contribution is simply \(w(x_i)x_i\psi(e_i)\).

A bounded \(\psi\) limits the residual part of a contribution. It does not limit the product with \(x_i\). Thus an observation can have an arbitrarily large effect on the estimating equation as its predictors grow. A leverage weight can control this product; one sufficient design is to choose \(w\) so that \(w(x)\|x\|\) is bounded. Such weighting must be designed with care: it changes the estimating equation and hence its population target, and it does not automatically make every form of contamination harmless.

Worked Example: A Response Outlier at an Ordinary Predictor Value

Consider a scalar-response regression with observations \((x,y)=(-1,0),(0,0),(0,10),(1,0)\). Treat the point \((0,10)\) as a response outlier. Use the Huber score \(\psi(r)=\max(-1.5,\min(r,1.5))\), a fixed scale of 1, and the estimating equations \(\sum_i\psi(r_i)=0\) and \(\sum_i x_i\psi(r_i)=0\), where \(r_i=y_i-a-bx_i\).

At \(a=0.5,\ b=0\), the three noncontaminated residuals are \(-0.5\), and the outlier residual is \(9.5\). Their scores are therefore \(-0.5,-0.5,-0.5,1.5\). The intercept equation is satisfied because \(-0.5-0.5-0.5+1.5=0\). The slope equation is also satisfied: the contributions from \(x=-1\) and \(x=1\) are \(0.5\) and \(-0.5\), while both points at \(x=0\) contribute zero.

For comparison, the least-squares intercept is \(\bar y=2.5\), and its slope is zero by symmetry. The robust fit has limited the outlier’s effect on the intercept. Its predictor value, however, is \(x=0\), so this example does not test resistance to high leverage.

Worked Example: A Bounded Residual Score With High Leverage

Use the scalar-response Huber score above in a regression with predictor vector \(x=(1,z)^T\), and examine the candidate line with intercept 0 and slope 1. Suppose a contaminating case has \(z=100\) and \(y=0\). Its residual is \(-100\), so its score is \(\psi(-100)=-1.5\).

Without predictor weighting, its estimating-equation contribution is \(\psi(-100)x=(-1.5,-150)^T\). Although the residual score has been clipped, the slope component has magnitude 150 because it is multiplied by \(z=100\).

Now use \(w(x)=\min(1,2/\|x\|)\). Since \(\|x\|=\sqrt{1+100^2}=\sqrt{10001}\), the weight is \(2/\sqrt{10001}\). The weighted contribution is \[ -\frac{3}{\sqrt{10001}}(1,100)^T. \] Its norm is exactly 3, since \((3/\sqrt{10001})\sqrt{10001}=3\). In contrast, the unweighted contribution has norm \(1.5\sqrt{10001}\), which is about 150. This bounds the case’s score contribution; it does not establish that the resulting fitted coefficient is unbiased or that its confidence region covers a clean-model coefficient.

What the Sandwich Covariance Establishes

The sandwich covariance describes the sampling variability of a regular estimating-equation root around its population target. It cannot turn a biased or inappropriate target into the desired coefficient. For a precise statement, let \(\theta=\operatorname{vec}(B)\), and write \(g(Z_i,\theta)\) for the vectorized score contribution, where \(Z_i=(Y_i,x_i)\). Suppose \(Z_1,\ldots,Z_n\) are independent and identically distributed, and define \(\theta_*\) by \(E[g(Z_i,\theta_*)]=0\).

Theorem (sandwich limit for an estimating-equation root): Suppose \(\theta_*\) is an interior, locally unique root; the estimator \(\hat\theta\) is consistent for \(\theta_*\); \(G=E[\partial g(Z_i,\theta_*)/\partial\theta^T]\) exists and is nonsingular; and \(E[g(Z_i,\theta_*)]=0\) with finite second moments. Assume the estimating function is differentiable near \(\theta_*\), and the usual uniform law of large numbers and local Taylor expansion hold there (in particular, the sample derivative converges to \(G\)). Then \[ \sqrt n(\hat\theta-\theta_*)\xrightarrow{d} N\!\left(0,G^{-1}\Omega G^{-T}\right),\qquad \Omega=E[g(Z_i,\theta_*)g(Z_i,\theta_*)^T]. \]

The usual sufficient regularity conditions for this conclusion include a central limit theorem for the score and a remainder in the sample Taylor expansion that is \(o_p(n^{-1/2})\). These conditions matter: the formula is not a distribution-free repair for arbitrary dependence, nonidentical contamination, or a nonregular estimator.

To derive the result, set \(P_ng(\theta)=n^{-1}\sum_{i=1}^n g(Z_i,\theta)\). Since \(P_ng(\hat\theta)=0\), a Taylor expansion around \(\theta_*\) gives

0=P_ng(\theta_*)+ \left\{P_n\frac{\partial g(Z_i,\tilde\theta)}{\partial\theta^T}\right\} (\hat\theta-\theta_*),

where \(\tilde\theta\) lies between \(\hat\theta\) and \(\theta_*\), with the expression understood through the integral form of the derivative if necessary. Consistency and the stated uniform convergence make the matrix in braces converge in probability to \(G\). The multivariate central limit theorem gives \(\sqrt n P_ng(\theta_*)\xrightarrow{d}N(0,\Omega)\), because the score has mean zero and finite second moments. Solving the expansion and applying Slutsky’s theorem yields

\sqrt n(\hat\theta-\theta_*)= -G^{-1}\frac{1}{\sqrt n}\sum_{i=1}^n g(Z_i,\theta_*)+o_p(1),

whose limiting covariance is \(G^{-1}\Omega G^{-T}\); the minus sign cancels in the covariance. A consistent plug-in estimate replaces \(G\) and \(\Omega\) with the sample derivative and the empirical outer product of the fitted scores. The covariance of \(\hat\theta\) is then estimated by \(n^{-1}\hat G^{-1}\hat\Omega\hat G^{-T}\).

Worked Example: Computing a Robust Sandwich Variance

Consider the intercept-only model \(Y_i=\beta+\varepsilon_i\), with known scale 1 and Huber score \(\psi(r)=\max(-1,\min(r,1))\). Suppose the errors independently take values \(-2,0,2\), each with probability \(1/3\). At the true coefficient \(\beta_*=0\), the three scores are \(-1,0,1\); their mean is zero, so \(E[\psi(Y_i-\beta_*)]=0\).

Except at the score’s two corner points, its derivative is 1 for \(|r|<1\) and 0 for \(|r|>1\). The error distribution assigns probability zero to those corners. Therefore \(G=E[-\psi'(Y_i-\beta_*)]=-1/3\) if the estimating function is \(g=\psi(Y_i-\beta)\). Equivalently, use \(A=-G=1/3\) in the positive-form sandwich. Also, \(\Omega=E[\psi^2]=2/3\). The asymptotic variance per observation is

G^{-1}\Omega G^{-T} =\frac{2/3}{(1/3)^2}=6.

Thus \(\operatorname{Var}(\hat\beta)\) is approximately \(6/n\), and the standard error is approximately \(\sqrt{6/n}\). In a sample containing exactly equal numbers of the three values, the empirical derivative magnitude is \(1/3\) and the empirical mean squared score is \(2/3\), giving the same sandwich calculation. This result quantifies variation around the estimating equation’s root; it does not show that a different contamination distribution would leave that root at \(\beta_*=0\).

Interpreting Inference Under Contamination

A key question is whether the population estimating equation remains centered at the coefficient being tested. Let \(B_0\) be the clean-model coefficient. If contamination changes the distribution so that \(E[g(Z_i,\operatorname{vec}(B_0))]\ne0\), then \(B_0\) is not the root \(\theta_*\) in the sandwich theorem. Under suitable regularity, the estimator may still converge to a root under the contaminated distribution, and the sandwich may estimate variability around that root. A confidence region calculated from it is not thereby a confidence region for \(B_0\).

The difference between the two targets is especially important for leverage contamination. A response score can be bounded while the predictor multiplier is not, so a few extreme predictor vectors can shift the estimating equation. If leverage weights are used, check what population equation those weights define and whether it is centered at the scientific coefficient of interest. Also check that the fitted derivative matrix is well-conditioned: a small eigenvalue can make the sandwich covariance large and coefficient contrasts unstable.

If the scale, residual scatter, or other nuisance quantities are estimated, treating them as fixed in the coefficient sandwich can omit their contribution to uncertainty. The standard approach is to write a joint estimating equation for the coefficient and nuisance parameters, then use its full derivative and score covariance. The covariance for the coefficient is obtained from the relevant block of that joint sandwich, not by silently substituting an estimated nuisance value into a formula that assumed it known.

Common Mistake and Diagnostic Questions

Common mistake: “The residual score is bounded, so no observation can dominate the robust regression, and the sandwich standard error accounts for any remaining contamination.” This reasoning overlooks the predictor multiplier: \(x_i\psi(e_i)\) need not be bounded. It also confuses variance estimation with protection against a shifted target. Check the full score contribution, the population estimating equation at the coefficient of interest, and the assumptions supporting the sandwich limit.

A practical diagnosis separates three questions. First, are unusual cases response outliers, high-leverage cases, or both? Second, does the estimating procedure bound or otherwise control their full score contributions? Third, does the covariance calculation estimate uncertainty for the coefficient being reported, under a sampling model that is credible for the data? A positive answer to only the second or third question is not enough to establish valid inference.

Check Your Understanding

Use the definitions and derivations above to answer the following questions.

Why can a bounded Huber residual score still produce an unbounded slope contribution in a regression with an extreme predictor?
For a leverage weight satisfying \(w(x)\|x\|\le C\) and a residual score satisfying \(\|\psi(e)\|\le c\), what bound follows for the norm of the weighted score contribution?
What population quantity does the sandwich covariance describe when the data follow a contaminated distribution and the estimating equation has a unique root?
Why does a sandwich standard error not, by itself, establish coverage for the clean-model coefficient?
If scale is estimated jointly with the regression coefficients, what change is needed in the sandwich calculation?