Residuals made independent

Residuals that look too normal found that a quantile plot of a regression’s residuals sees the errors’ shape less well the more columns the fit takes. Each residual is a mixture of every row’s error, mixtures of independent draws move towards the normal, and the straightness test that flags a skewed error in 82.9% of samples of forty errors flags it in 51.2% of the residuals of a ten-column fit and 23.0% of a twenty-column one. The loss was blamed on the mixing, and the essay ended on the obvious repair: residuals built so that each depends on one new error and a finite past, the recursive residuals, which are exactly independent and identically distributed under normal errors.

They are exactly that, and it does not help. A second construction that is optimal by design does not help either, and the reason is a fact about all such constructions at once.

Two ways to make residuals independent #

Recursive residuals take the rows in some order. For each row after the first $p$, fit the model to the rows before it, predict the row, and divide the prediction error by its standard deviation:

$$w_t = \frac{y_t - x_t^{\top} b_{t-1}}{\sqrt{1 + x_t^{\top}(X_{t-1}^{\top}X_{t-1})^{-1}x_t}}.$$

There are $n - p$ of them. Each is a fixed linear combination of the errors, and the combinations are orthonormal and orthogonal to the design’s columns, so under independent normal errors of one spread the $w_t$ are exactly independent normals of that spread — for any design and any order of the rows. A quantile plot of them can be read against the band for $n - p$ independent normals with no approximation at all.

BLUS residuals, Theil’s construction, choose $n - p$ of the rows to stand for and ask, among all sets of residuals with the same exact null distribution, for the one closest in expected squared distance to those rows’ own errors. With $M = I - H$ the ordinary residual-maker and $J$ the rows kept — here every row but the $p$ of highest leverage — the answer is $(JMJ^{\top})^{-1/2}JM$. No residuals that are independent under the null are nearer the errors they stand for.

Both keep their level exactly. Across five designs — a line with an evenly spread covariate, a line with one point at leverage 0.73, and fits of five, ten and twenty columns — the five shape tests on recursive residuals reject between 3.9% and 6.4% of samples with normal errors, which is sampling error around 5% on three thousand samples, and in every order of the rows tried.

What they see #

The ordinary residuals of a line flag a skewed error in 79.1% of samples; BLUS residuals in 76.4%; recursive residuals in 64.4%. At five columns, 71.4%, 64.9% and 49.5%. At ten, 51.2%, 38.3% and 25.9%. At twenty, 23.0%, 16.2% and 9.1%. At every design the order is the same: the ordinary residuals, with all their mixing and all their correlation, see the shape best, the optimal independent set second, and the recursive residuals worst by a wide margin.

The other departures agree. Against light tails read by the kurtosis test, ten-column residuals flag 17.7% of samples, BLUS 12.1% and recursive 8.0%; against heavy tails read by the worst-point test, twenty-column residuals flag 13.2%, BLUS 10.2% and recursive 7.1%. The one near-exception is light tails at twenty columns, where everything sits near the level — 4.7%, 6.3% and 5.0% — because nothing there sees anything.

The benchmark the repair was meant to approach is the dashed line: $40 - p$ independent errors, what the fit would cost if it only took observations away. At ten columns that is 69.1%; at twenty, 49.5%. Neither independent construction comes near it.

Every independent set is a rotation #

The reason is that the independent residuals were never going to contain anything the ordinary ones do not. Write the recursive residuals as $w = C\varepsilon$. Because the rows of $C$ are orthonormal and span exactly the space the ordinary residuals live in — everything orthogonal to the design’s columns — the residual-maker factors as

$$M = C^{\top}C,$$

so the ordinary residuals are $M\varepsilon = C^{\top}w$: the recursive residuals, rotated back. The same holds for BLUS residuals and for every other set with an exact independent null. They are all the same $n - p$ numbers’ worth of information about the errors, expressed in different coordinates, and a fit that took $p$ dimensions of the errors away took them from all of them equally.

What differs between the coordinate systems is how the information about each error is spread across the residuals, and a shape test is a nonlinear reading of the residuals’ marginal distribution, which cares a great deal about that spread. The ordinary residuals keep as much of each error in its own residual as the fit allows. The independent constructions have to redistribute it.

Independence is bought by mixing #

For an ordinary residual the share of its variance that is its own error is $1 - h_{ii}$, and it averages $(n-p)/n$ exactly: 0.750 at ten columns on forty rows, 0.500 at twenty. BLUS residuals average 0.711 and 0.498 — nearly as pure as the ordinary ones, but only over the $n - p$ rows kept, so the $p$ rows left out contribute to every one of them and to none as their own. Recursive residuals average 0.562 and 0.294.

The recursive average hides its shape. A recursive residual early in the sequence predicts its row from barely more rows than there are columns, the prediction is poor, and the residual is mostly the fitted rows’ errors: the least pure recursive residual at twenty columns is 0.4% its own error, and at ten, 7.3%. Late in the sequence the prediction is good and the residual is nearly pure. Independence across residuals is achieved exactly by piling the mixing into the early ones, and a mixture is what a shape test cannot see through, which is the whole of the first essay’s finding.

So the question the last essay asked has a precise answer. Mixing costs power, but correlation between residuals was never the mechanism: residuals can be made uncorrelated only by mixing some of them harder, and the ordinary residuals — correlated, unequal, and as unmixed as the fit permits — are the best of the residual sets for reading a shape.

One sample, read four ways #

The two effects can be seen in a single sample. Take forty lognormal errors and a twenty-column fit — the sample the first essay used to show a skew being absorbed. The errors themselves have a straightness statistic of 0.1259 and a skewness of 1.161: plainly skewed, well past the 5% points for forty observations of 0.0551 and 0.680. The ordinary residuals of the fit read 0.0259 and 0.276, well inside them.

The recursive residuals of the same errors read 0.0461 and 0.440, and the BLUS residuals 0.0385 and 0.529. Each is further from normal than the ordinary residuals — the independent constructions have not hidden the skew as thoroughly — and neither crosses its line, because there are only twenty of them and the 5% points for twenty observations are 0.0963 and 0.866. An independent set is read against the band for the number of values it has, and it has $n - p$.

So the independent sets pay twice. They carry the information of $n - p$ values and are read as $n - p$ values, with the wider band that implies; and they are more mixed on average than the ordinary residuals, which blunts even those values. The ordinary residuals carry the same $n - p$ dimensions of information spread over forty values, and — as the first essay found — their forty values are close enough to forty independent normals under the null that the band for forty holds its level. Reading a rotation of twenty dimensions as forty values is what an exact construction gives up, and at these sizes it is worth more than exactness.

The comparison with twenty pure errors separates the two costs. Twenty independent lognormal errors are flagged by the straightness test 49.5% of the time with the band for twenty; the twenty recursive residuals of a twenty-column fit, with the same band, 9.1%. Everything between those two numbers is mixing, and it is most of the loss.

The order a recursion takes #

Recursive residuals depend on the order of the rows, and for a design with a far point the dependence is large.

On a line with one point at leverage 0.73, recursive residuals taken with the far point first see a skewed error 75.1% of the time; with it last, 50.1%. First, the far point is among the rows the recursion starts from and never gets a residual of its own; last, its residual is a prediction far outside the range of every row before it, and that prediction’s error is mostly the fitted line’s errors carried out to the far point. At ten columns the spread is smaller — 31.2% with high leverage first, 18.6% with it last, 25.9% and 24.9% in the natural and a shuffled order — but it is there in every design. The level holds in every order; the power is a property of the order chosen.

That is a strange property for a diagnostic to have, since nothing about the errors’ shape depends on which row was collected first, and it is the practical reason recursive residuals are not used to read a shape. What they are used for is the thing their order-dependence suits: a test for a change in the relation partway through an ordered sample, where a run of large recursive residuals after some row says the rows after it are not predicted by the rows before. That is a question about the order, and a break that was looked for is where that question has been measured.

Where the recursion puts the information #

The recursive construction has one more property worth stating, because it is the reason the order matters so much. The information the fit leaves about the errors is spread by the recursion unevenly along the sequence: a late residual is nearly its own row’s error, an early one is nearly all the fitted rows’ errors. Any feature of the errors that lives in a particular row — a skewed error that happened to fall on the far point, one wild value — is seen clearly if that row comes late and hardly at all if it comes early.

That is why, with the far point last, the recursion does worse than in any other order: the far point’s own residual is dominated by the fitted line carried out to it, so the row whose error most affects the fit contributes least of its own error to the plot. With the far point first it is among the rows the recursion starts from and has no residual at all, and the other thirty-eight residuals are spread over rows whose predictions are ordinary. Neither choice is wrong in the level; the second is simply the order in which the most information arrives late.

A reader cannot know which order a software package chose, and different packages choose differently — the order of the data file, the order of time, the order of a covariate. A shape diagnostic whose power depends on a choice the reader cannot see is weaker than its exact level suggests, and twenty residual plots already found that a reader has little idea what a correct plot of forty ordinary residuals looks like, let alone one whose construction changes with the rows’ order.

What reading a residual plot comes to #

Residuals are not the errors showed that the ordinary residuals’ spreads differ by the design; the band the eye was standing in for and a number for the shape showed what a quantile plot’s band and its one-number summaries can and cannot see on forty independent observations. Put together with the finding here, the reading of a residual plot comes to four statements.

Read the ordinary residuals. Of every set of residuals that has an exact or near-exact null, they see an error’s shape best, because they keep as much of each error in its own residual as the fit permits.

Read them against the band for n, or a band simulated from the design. The first essay found the two within a point of each other; the recursive residuals’ exact band buys exactness in a construction that sees less.

Count the columns before believing a clean plot. At twenty columns on forty rows, no construction of residuals sees a skewed error more than a quarter of the time, because the fit took half of the errors’ dimensions and no rotation of what is left restores them. $R^2$ is not a fit found that twenty useless predictors on thirty points manufacture an $R^2$ of 0.69 from noise; the same columns manufacture a normal-looking residual plot from skewed errors.

State the share of the errors the fit leaves before the plot is drawn. The ordinary residuals’ own-error share averages $(n-p)/n$, and it is known from the design alone, before any outcome is measured, exactly as the points a slope rests on found the effective number of points behind a slope to be. A plot drawn from residuals that are half their own errors on average is a plot of a half-strength signal, and a reader told that number in advance knows how much a clean plot can mean. At two columns on forty rows the share is 0.95 and the plot is nearly a plot of the errors; at twenty it is 0.50, and no reading of the residuals, rotated or not, will restore the other half.

Use recursive residuals for the question they answer. They are the right residuals when the question is whether the relation changed along the order of the rows, and the wrong ones for whether the errors are normal.

Exact, and counted #

Exact: recursive and BLUS residuals of normal errors are independent normals for any design and any order, and every set of residuals with that property is a rotation of the ordinary residuals, $M = C^{\top}C$. The ordinary residuals’ own-error share averages $(n-p)/n$.

Counted, over three thousand samples of forty with lognormal errors of σ = 0.5: the straightness test on ordinary, BLUS and recursive residuals flags 51.2%, 38.3% and 25.9% of samples at ten columns and 23.0%, 16.2% and 9.1% at twenty; recursive residuals of a line with one far point flag 75.1% with the far point first and 50.1% with it last.

Not claimed: that no residual construction can do better than the ordinary one. A nonlinear construction — residuals from a robust fit, or residuals rescaled row by row by $1/\sqrt{1 - h_{ii}}$ — escapes the rotation argument, which covers only linear sets with an exact independent null. The studentised residuals many packages plot are one such rescaling; whether any nonlinear construction reliably beats the ordinary residuals on shape is not measured here.

Still open: residuals of a model that was chosen #

Every fit here has its columns fixed before the data arrive. In practice the columns are often chosen from the data — a stepwise search, a lasso, a set of interactions kept because they were significant — and a column chosen because it fits the errors takes more of the errors’ shape with it than a column fixed in advance. A variable that happens to line up with the three largest errors is the variable a search keeps, and the residuals it leaves have had their right tail removed on purpose.

How much more normal the residuals of a selected model look than those of a fixed model of the same size, whether a band simulated by repeating the selection on simulated errors restores the level, and whether a reader told only the final model’s column count is reading a plot whose power is that of a much larger fit, are the measurements the residual plots of most published regressions would need. None of them has been made here.