Why Residual Exchangeability Needs to Be Checked
Residual permutation can provide a reference distribution for a test involving a predictor while nuisance covariates remain in the model. Its appeal is practical: fit the nuisance model, permute its residuals, and use the permuted data to represent outcomes under a null hypothesis. But a permutation test is justified only when the objects being permuted have the required invariance. In particular, residuals that differ systematically in variance or are strongly affected by leverage need not be exchangeable.
Here, exchangeability means that the joint distribution does not change when the components are permuted. It is stronger than saying that all residuals have mean zero, and stronger than saying that their average variances are similar. In genomic analyses, a common design matrix may make leverage patterns repeat across many features, while the error variances may also vary by sample, group, or fitted value. Those patterns can make a residual-permutation test unreliable even when the fitted mean model appears reasonable.
The earlier material on residual permutation with nuisance covariates motivates the procedure. This tutorial focuses on a diagnostic question: do the residuals have a distribution for which the planned permutations are credible?
Residual Variance, Leverage, and Heteroskedasticity
Consider the linear model \(y=X\beta+\varepsilon\), where \(y\) is an \(n\)-vector, \(X\) is a fixed design matrix of full column rank, and \(\beta\) is a vector of regression coefficients. Assume \(E(\varepsilon\mid X)=0\). The ordinary least squares (OLS) fitted values and residuals are \(\hat y=Hy\) and \(e=y-\hat y=My\), where \(H=X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\) is the hat matrix and \(M=I-H\). The diagonal entry \(h_{ii}\) is the leverage of observation \(i\).
If errors are homoskedastic and uncorrelated, then \(\operatorname{Var}(\varepsilon\mid X)=\sigma^2 I\), so \(\operatorname{Var}(e\mid X)=\sigma^2M\). In particular, \(\operatorname{Var}(e_i\mid X)=\sigma^2(1-h_{ii})\). Thus even with a common error variance, residuals at high-leverage observations have smaller marginal variance. Also, residuals are generally correlated because the off-diagonal elements of \(M\) need not be zero.
Under heteroskedastic, uncorrelated errors, write \(\operatorname{Var}(\varepsilon\mid X)=\Omega=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2)\). Then \(\operatorname{Var}(e\mid X)=M\Omega M^\mathsf{T}\), and the individual residual variances are \(\operatorname{Var}(e_i\mid X)=\sum_{j=1}^n M_{ij}^2\sigma_j^2\). The simple homoskedastic expression \(\sigma^2(1-h_{ii})\) no longer applies. Both unequal error scales and the design-dependent residual transformation affect the distribution of the residual vector.
A useful sufficient condition for exchangeability of a mean-zero multivariate normal residual vector is that its covariance matrix be invariant under the intended permutations: \(P\operatorname{Var}(e\mid X)P^\mathsf{T}=\operatorname{Var}(e\mid X)\). For a normal vector, the mean and covariance determine its joint distribution. With nonnormal errors, covariance invariance alone is not enough; the full joint distribution must be invariant. Neither condition follows automatically from fitting an OLS model. Consequently, these calculations diagnose important failures, but passing a marginal-variance check does not prove that a residual permutation test is exact.
Worked Examples: Calculating Leverage and Residual Compression
Worked Example: A High-Leverage Observation
Suppose a simple regression with an intercept has five predictor values \(x=(8,0,1,2,3)\). Their mean is \(\bar{x}=2.8\), and the sum of squared deviations is \(S_{xx}=38.8\). Under homoskedastic, uncorrelated errors with variance \(\sigma^2\), calculate the leverage and residual variance for the observation with \(x_i=8\).
For simple regression with an intercept, \(h_{ii}=1/n+(x_i-\bar{x})^2/S_{xx}\). Therefore,
The arithmetic matters: \(27.04/38.80\approx0.6969\), so the leverage is approximately \(0.8969\), not \(0.8979\). The residual variance is
If \(\sigma=2\) response units, the residual standard deviation at this design point is approximately \(2\sqrt{0.1031}=0.6422\) response units. The small residual scale does not show that this observation has a small error: high leverage has made its fitted value strongly dependent on its own response, leaving a compressed residual. The calculation assumes the stated homoskedastic, uncorrelated error model; it does not establish that residuals from this design are exchangeable.
Worked Example: Unequal Error Variances Change Residual Scales
Worked Example: Residual Variances in an Intercept-Only Fit
Consider three independent, mean-zero Gaussian errors with variances \(1,1,\) and \(9\) response-units squared. Fit an intercept only. Calculate the variance of each OLS residual and assess whether the residuals are exchangeable.
For an intercept-only fit, \(H\) has every entry equal to \(1/3\), so \(M=I-H\) has diagonal entries \(2/3\) and off-diagonal entries \(-1/3\). The residual for observation 1 is \(e_1=(2\varepsilon_1-\varepsilon_2-\varepsilon_3)/3\). Independence gives
The same calculation gives \(\operatorname{Var}(e_2\mid X)=14/9\). For observation 3, the coefficient on its own error is \(2/3\), while the other two coefficients are \(-1/3\), so
The third residual has a larger marginal variance than the first two, so the residual vector cannot be invariant to all permutations: swapping the first and third components would swap unequal marginal variances. The residuals also sum to zero because the intercept was fitted, which creates dependence. This example shows why “the model included an intercept” and “the residuals average to zero” are not evidence of exchangeability.
Worked Example: Diagnosing a Variance Pattern
Worked Example: Residual Spread Across Fitted Values
Suppose a hypothetical feature is measured in six samples. After fitting the null model that includes the nuisance covariate, the fitted values are \(2,2,2,8,8,8\) expression units. The corresponding residuals are \(-1,0,1,-4,0,4\) expression units. Compare the residual spread across the two fitted-value levels.
At fitted value 2, the residual mean is \((-1+0+1)/3=0\), and the sum of squared deviations is \(1+0+1=2\). The sample variance is \(2/(3-1)=1\), in squared expression units. At fitted value 8, the residual mean is also zero, and the sum of squares is \(16+0+16=32\); the sample variance is \(32/2=16\) squared expression units. The second group's sample standard deviation is therefore \(4\) expression units, compared with \(1\) expression unit in the first group.
This fourfold standard-deviation difference is a clear warning that residuals from the two fitted-value regions may not be interchangeable. The groups are small, and the residuals are themselves fitted quantities, so these sample variances are diagnostic evidence rather than proof about the unknown error variances. The pattern calls for checking the measurement process and the model, not for treating the apparent variance ratio as a definitive estimate.
A Practical Diagnostic Sequence
Diagnostics should be performed on the exact model and set of observations used to generate the residuals that will be permuted. In a genomic scan, the design-related leverage can be common across features, but residual variance patterns can differ by feature. A single reassuring plot or an overall residual standard deviation does not certify exchangeability.
Calculate \(h_{ii}\), inspect its distribution, and identify unusually influential design points. In an intercept-only model all leverages are \(1/n\); with covariates, leverage can be concentrated in a few observations. A high leverage value is a warning about residual compression, not a standalone proof of invalid inference.
Plot residuals against fitted values and important covariates, and compare spread across groups or batches. A funnel shape or systematically different group spreads is evidence against a common error scale. Check that a few extreme residuals are not driving the apparent pattern.
Under homoskedastic errors, \(e_i/\sqrt{1-h_{ii}}\) has variance \(\sigma^2\), when \(h_{ii}<1\). This can reveal leverage-related scale differences, but does not remove residual dependence or prove joint invariance. Under heteroskedasticity, this adjustment generally does not equalize variances.
Ask which observations the procedure exchanges and whether their error distributions are plausibly alike under the null. If variance or leverage patterns distinguish observations that the algorithm freely swaps, report the concern and use a method with calibration appropriate to the design and error structure.
What the Diagnostics Can and Cannot Establish
The residual vector is produced by the projection \(M\), so equal variance of the original errors would not make the residuals independent. High leverage creates particularly unequal residual variances under homoskedasticity; heteroskedasticity adds unequal error scales before projection. When the resulting joint distribution is not invariant under the proposed permutations, the permutation reference distribution need not match the null distribution of the test statistic. The resulting p-values may therefore be miscalibrated; the direction and size of the error are not determined by a diagnostic plot alone.
A common correction to the reasoning is to standardize residuals by their estimated standard deviations and then assume they can be permuted freely. Standardization may address a marginal scale difference under a correctly specified variance model, but it does not generally remove dependence or other differences in joint distribution. Estimated standardizers also add uncertainty. Such transformations require their own justification and are not a general guarantee of an exact permutation test.
Residual permutation is also different from randomization inference based on an experiment's actual assignment mechanism. If treatment labels were randomized, that assignment mechanism can justify permutations of labels under a sharp null, subject to the design restrictions. It does not automatically justify permuting fitted residuals. Model-based alternatives include inference using heteroskedasticity-robust standard errors, with their assumptions and finite-sample limitations assessed for the setting. Wild bootstrap methods are another approach for nonexchangeable errors; they require their own calibration choices and are not made valid merely by observing unequal variances.
In high-dimensional testing, marginal calibration matters before any multiple-testing adjustment is applied. A correction for multiplicity cannot generally repair p-values that are systematically invalid because the null reference distribution is wrong. Conversely, a flagged diagnostic does not prove that every resulting p-value is wrong. The appropriate conclusion is that residual permutation needs stronger justification or a more suitable calibration strategy.
Common Mistakes
- Equating zero residual mean with exchangeability. OLS residuals satisfy constraints such as \(X^\mathsf{T}e=0\), but those constraints do not make their joint distribution permutation-invariant.
- Using the wrong leverage arithmetic. In the five-point example, \(27.04/38.80\approx0.6969\), yielding \(h_{ii}\approx0.8969\) and residual variance \(0.1031\sigma^2\). Rounding an intermediate calculation incorrectly changes every subsequent value.
- Assuming residual variance equals error variance. Even under homoskedasticity, \(\operatorname{Var}(e_i\mid X)=\sigma^2(1-h_{ii})\), not \(\sigma^2\).
- Treating a diagnostic plot as a formal proof. A plot can reveal serious warning signs, but small samples have limited power to reveal variance patterns, and a visually uniform spread does not establish joint exchangeability.
- Assuming multiplicity adjustment repairs miscalibration. Multiple-testing procedures cannot be relied on to restore valid individual null p-values when the permutation reference distribution is inappropriate.
Check Your Understanding
Answer each question, then compare with the concise solutions.
- Under homoskedastic, uncorrelated errors, what is \(\operatorname{Var}(e_i\mid X)\) in terms of \(h_{ii}\) and \(\sigma^2\)? What happens as leverage approaches 1?
- Why do unequal marginal residual variances rule out invariance under a permutation that swaps the corresponding residuals?
- For an intercept-only model with three independent errors of variances \(4,4,4\), what are the three residual variances?
- Does dividing residual \(e_i\) by \(\sqrt{1-h_{ii}}\) prove that residuals can be freely permuted? Explain.
- Name one diagnostic pattern that would make residual permutation concerning, and state what the pattern does and does not establish.
Solutions. (1) The variance is \(\sigma^2(1-h_{ii})\); it approaches zero as \(h_{ii}\) approaches 1, reflecting residual compression. (2) A permutation-invariant joint distribution must have the same marginal distribution at positions that can be swapped, including the same variance; unequal marginal variances contradict that requirement. (3) In this homoskedastic case each leverage is \(1/3\), so each residual variance is \(4(1-1/3)=8/3\) squared response units. (4) No. Under the stated homoskedastic model the division equalizes marginal variances, but residual dependence and joint distribution invariance remain to be justified; with heteroskedastic errors the adjustment may not even equalize marginal variances. (5) A funnel-shaped residual-versus-fitted plot is one warning sign. It suggests nonconstant residual spread and questions the common-scale assumption, but by itself does not quantify the null p-value error or prove exactly how the test is miscalibrated.