Why don't standard estimation methods consider the distribution of the residual?

Question

Consider a linear model in which a vector of data $\mathbf{d} \in \mathbb{R}^M$ is related to an unknown parameter of interest $\mathbf{x} \in \mathbb{R}^N$ via $$\mathbf{d} = \mathbf{A}\mathbf{x} + \mathbf{n},$$ where $\mathbf{A} \in \mathbb{R}^{M\times N}$ is a known matrix with rank $N$ and $\mathbf{n} \in \mathbb{R}^M$ is an unknown vector of zero-mean i.i.d. Gaussian noise realizations with known variance $\sigma^2$ (i.e., $\mathbf{n} \sim \mathcal{N}(\mathbf{0},\sigma^2\mathbf{I})$).

It is well known that the maximum likelihood estimator of of $\mathbf{x}$ is also the least-squares estimator and (via the Gauss-Markov theorem) the best linear unbiased estimator (BLUE): $$\hat{\mathbf{x}}_{\mathrm{ML}} = \arg\min_\mathbf{x} \|\mathbf{A}\mathbf{x}-\mathbf{d}\|_2^2= (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{d}.$$

To me, it seems intuitive to interpret the residual of this maximum likelihood estimator $$\mathbf{r}_{\mathrm{ML}} = \mathbf{d} - \mathbf{A}\hat{\mathbf{x}}$$ as some form of estimate of the Gaussian noise $\mathbf{n}$.

However, there are cases where the empirical histogram of $\mathbf{r}_{\mathrm{ML}}$ will be wildly different from the expected histogram of a vector $\mathbf{n}$ generated by a zero-mean i.i.d Gaussian distribution with known variance $\sigma^2$. For example, if we let $\mathbf{A}$ be the identity matrix, then $\hat{\mathbf{x}}_{\mathrm{ML}} = \mathbf{d}$, and $\hat{\mathbf{r}}_{\mathrm{ML}} = \mathbf{0}$. The empirical histogram of $\hat{\mathbf{r}}_{\mathrm{ML}}$ has a point mass at 0, which is very different from the empirical histogram one would expect from i.i.d. Gaussian noise.

This is nonintuitive to me, and I have two questions:

1. Is it reasonable to want the empirical histogram of the residual vector $\hat{\mathbf{r}}_{\mathrm{ML}}$ to match the distributional prior information I have about the noise $\mathbf{n}$?
2. Why don't standard estimators (like the maximum likelihood estimator) consider the empirical histogram of the residual vector? It seems like it would be easy to construct a penalized maximum likelihood method that penalizes deviations between the residual and the expected empirical histogram.

Answer 1

The first thing to ask is what the aim is of estimating a regression hyperplane by Least Squares (which is Maximum Likelihood for a model with normal errors). The aim is often to fit the data as well as possible. When doing model-based inference, the aim is to recover a hyperplane that "explains" the responses (here called ${\bf d}$) best in terms of the explanatory variables (here stored in ${\bf A}$). Normality is a model assumption for the underlying error term. Note that normality is required for an optimality statement and precise distributions of certain standard tests used in regression, but it is not required for the regression to make sense (in the sense of fitting the data well; even distribution theory and tests will in very many non-normal situations be approximately valid). The distribution of residuals is often used to diagnose normality of the error term, even though it won't be perfectly normal even if the error term is, see below. A basic condition of (frequentist) model-based statistics is that we cannot directly observe whether model assumptions are fulfilled (arguably they never are; "all models are wrong but some are useful"). In particular, there is no assumption that the empirical distribution of residuals is normal, and defining a "penalty" as mentioned in the question for enforcing it to look more normal does not help with actually fulfilling the normality assumption in the model. In terms of diagnostics, the residuals are more useful if they are not forced to look as normal as possible (which of course would have an effect even if the underlying true distribution were not normal).

I have added the introduction above after some time because I thought I should say more about the "role" of normality in regression, and how the idea in the question doesn't help matters. Here is my original response:

The standard situation for a linear model is that the number of explanatory variables or rather, the corresponding regression parameters ($N$ in your notation) is smaller than the number of observations ($M$ in your notation, which is generally nonstandard; I just mention this for readers who get confused by the choice of notation). If ${\bf A}$ is the identity matrix, $N=M$, and indeed the observations can all be fitted perfectly (residual zero) by the regression. Note that there is the concept of "degrees of freedom" of the residuals. The residuals have $M-N$ degrees of freedom, meaning that the larger $N$ is compared to $M$, the stronger is the empirical residual distribution determined by the parameter estimation. For example, with $M=1$, the empirical mean of residuals must be exactly zero (note that this is not normally the case for $M$ normally distributed random variables with theoretical mean zero, as these will have a random empirical mean usually close to but not exactly equal to zero); if $M=N$ all residuals will be zero (perfect fit), assuming that ${\bf A}$ is of full rank.

Note that this even holds if indeed ${\bf n}$ is perfectly normally distributed as assumed. So the empirical residuals $\hat{\bf r}_{ML}$ will in this way not perfectly match the distribution of ${\bf n}$.

For estimating ${\bf x}$, we want the observations to be fitted well in the first place, i.e., we want the residuals to be generally small (in absolute value, or squared). In this sense, $\hat{\bf r}_{ML}=0$ is as good as it gets; we wouldn't improve matters if we tried to make $\hat{\bf r}_{ML}$ looking more normal.

In fact, we may also use Least Squares estimation of ${\bf x}$ in case that ${\bf n}$ is not normally distributed but has mean 0. This is maximum likelihood for normal ${\bf n}$, but is also reasonable for other distributions for which it is not maximum likelihood, as the Least Squares criterion measures the quality of approximation of the observations by the estimated regression in a rather general way. This is another reason why we wouldn't want to enforce the residuals to look more normal at the cost of losing fit quality; in reality we can't know the true distribution of ${\bf n}$, which may well be non-normal (this is the reason why we can't just say, if ${\bf n}$ has a non-normal distribution, we take the maximum likelihood estimator for that distribution), but Least Squares may still give us a good regression hyperplane. So your idea of penalising the estimator so that the residual histogram looks more normal will not improve the regression estimator regarding the thing most users want to achieve in the first place.

You may be aware that in fact penalisation is often used in regression if $N$ is large compared to $M$, including $M=N$ as in your example ("Lasso"). The aim of this penalisation, however, is not to make the residuals look more normal, but rather (motivated by theory of overfitting) to shrink the absolute values of the components of $\hat {\bf x}$ and reduce some of them to zero (variable selection).

Answer 2

1. Is it reasonable to want the empirical histogram of the residual vector $\hat{\mathbf{r}}_{\mathrm{ML}}$ to match the distributional prior information I have about the noise $\mathbf{n}$?

Yes that is unreasonable. That is why instead we consider the sampling distribution of the residuals to be differently distributed from the errors.

$$\mathbf{r} = \left(\mathbf{I}-\mathbf{A}(\mathbf{A}^t\mathbf{A})^{-1}\mathbf{A}^t\right) \mathbf{y}$$

In the case that it makes sense to describe the distribution of the outcome as a sum with zero mean error terms (such as with normal distributed errors) and a bias term $$\mathbf{y} = E(\mathbf{\hat{y}}) + {\mathbf{y}}_{\text{bias}} + \mathbf{n}$$

then

$$\mathbf{r} = \left(\mathbf{I}-\mathbf{A}(\mathbf{A}^t\mathbf{A})^{-1}\mathbf{A}^t\right) \mathbf{n} + \left(\mathbf{I}-\mathbf{A}(\mathbf{A}^t\mathbf{A})^{-1}\mathbf{A}^t\right) {\mathbf{y}}_{\text{bias}} $$

or if there is no bias (when the true model is in the matrix $\mathbf{A}$) then

$$\mathbf{r} = \left(\mathbf{I}-\mathbf{A}(\mathbf{A}^t\mathbf{A})^{-1}\mathbf{A}^t\right) \mathbf{n}$$

An example for the way that people incorporate this is by estimating the variance of the noise $\mathbf{n}$ by considering the sum of squared residuals divided by the number $n-p$ instead of $n$ (where $n$ is the number of observations, and $p$ the rank of the matrix, the number of parameters). In your case with the identity matrix as design matrix, this will lead to the extreme case that the residuals will have zero variance.

2. Why don't standard estimators (like the maximum likelihood estimator) consider the empirical histogram of the residual vector?

You also have estimators that use the residuals in some way such as some robust estimators which fit multiple times and adjust weights based on the residuals of earlier fits, or methods that estimate the median instead of the mean. Also, in the analysis of the statistical properties of the estimator, people consider the residuals when they apply methods like Monte Carlo simulations.

It seems like it would be easy to construct a penalized maximum likelihood method that penalizes deviations between the residual and the expected empirical histogram.

The point/principle of this penalty is not so much explained. What would be the effect of such penalty? Maximumizing the likelihood is already the method to minimize the deviations between the empirical distribution and the assumed distribution.

In a way you can see the maximum likelihood as minimizing the cross entropy between the empirical distribution and the expected distribution.

Consider an empirical density/mass distribution that is the derivative of the empirical distribution $\hat{f}(x) = \frac{\text{d}}{\text{d}x} \hat{F}(x)$, by considering the derivative of the step function as Dirac delta functions

$$\hat{f}(x) = \frac{1}{n} \sum_{i=1}^n \delta_{r_i}(x)$$

then the maximum likelihood method is equivalent to minimizing the cross-entropy between empirical $\hat{f}(x)$ and expected $f(x)$

$$\begin{array}{rcl} H(\hat{f},f) &=&- \int \hat{f}(x) \log f(x) \text{d}x \\ &=&- \int \frac{1}{n} \sum_{i=1}^n \delta_{r_i}(x) \log f(x) \text{d}x \\ &=& -\frac{1}{n} \sum_{i=1}^n \delta_{r_i}(x) \log f(x) \text{d}x \\ &=& -\frac{1}{n} \sum_{i=1}^n \log f(r_i)\\ &\propto&- \log \prod_{i=1}^n f(r_i) \end{array}$$