How to select the "best" distribution of the errors in linear data?

Question

The maximum likelihood estimator for linear data $y_i = \beta x_i + \epsilon_i$ with i.i.d. normally distributed errors $\epsilon_i \sim \mathcal{N}(0,\sigma^2)$ is ordinary least squares, i.e.

$$\hat{\beta}_{\text{Normal}} = \text{argmin}_\beta \sum_i (y_i-\beta x_i)^2.$$

However, if the errors are Laplace distributed, $\epsilon_i \sim \mathcal{L}(0, b)$, the MLE turns out to be

$$\hat{\beta}_{\text{Laplace}} = \text{argmin}_\beta \sum_i |y_i-\beta x_i|,$$

which is L1 regression. In principle, I don't see anything preventing our errors from following any zero-mean distribution, and deriving an appropriate minimization problem, either analytically or computationally.

However, how can we find the distribution of the errors? For example, suppose we fit an L1 regression to the data and a statistical test confirms the residuals $r_i$ are Laplace distributed. Then we can be reasonably confident that we selected the "right type of regression" (L1 regression rather than, say, OLS). But what if we perform an L1 regression and the residuals $r_i$ don't look Laplace distributed?

In general, is there some sort of procedure for selecting the "best" distribution for $\epsilon_i$, in the sense that when you solve the corresponding minimization problem for $\hat{\beta}$, the residuals $r_i = y_i - \hat{\beta}x_i$ follow the same distribution you initially assumed the errors $\epsilon_i$ would follow?

Answer 1

You can always try different distributions and compare (penalized) likelihoods to get a measure of what fits best. This provides a solution, but can be criticized for being purely data-based as opposed to theory-based. Also many models will fit nearly equally well, so it is not at all definitive. You also have to take care when comparing models with transformed Y versus non-transformed Y, and you have to take care when comparing discrete versus continuous models. Nevertheless, it provides an answer, and there are several worked examples in here https://www.google.com/books/edition/Understanding_Regression_Analysis/y3jtDwAAQBAJ?hl=en using this approach.

One bad idea is to fit a distribution to the residuals. While you will find many people on this site stating that the distribution assumption refers to the errors, and not to Y, this statement can be misleading and lead to bad practice. It is misleading in the sense that the assumption really refers to the conditional distributions of the errors, which equivalently means that it refers to the conditional distributions of the Y's.

So if one were to fit a distribution to the residuals in the hypothetical case where the conditional distributions were normal but heteroscedastic, one would likely conclude that the distribution was heavy-tailed rather than normal.

Further, if the Y data were discrete, and a distribution was found using the residuals, a continuous model would be chosen incorrectly.

Answer 2

The first thing to understand about issues like this one is that statistical models live in the land of mathematics rather than in real life, and in real life no statistical model is ever fulfilled (arguably frequentist probabilities don't really exist in the first place, let alone i.i.d. residuals, continuous distributions, and whatnot).

So in any real situation, neither will the normal be the "true" error distribution, nor the Laplace, nor actually anything that can be written down.

The reason why we use model-based theory is that it gives us an idea what to do, which will work well in an idealised situation in the world of mathematics. This is is at least something, even though it doesn't guarantee us that the method will work well in reality. But the idea that in order to apply a method derived from a model-based approach, the model must be true in reality is wrong (if that were the case, we could never apply anything).

Both the LS-regression and the L1-regression (the ML for Laplace) fit a regression line/hyperplane to the data in a well defined and intuitively reasonable way (minimising the sum of squared, or unsquared deviations). When it comes to choosing between them (or even other methods such as proper breakdown robust estimators such as MM), we need to use all understanding that we have, which includes model-based theoretical results, but also other things, like regarding the way how outliers affect both of these estimators (the LS somewhat more), numerical instability of L1, potentially good asymptotic properties of the LS-estimator and derived inference for many non-normal distributions due to the central limit theorem.

Obviously we may hope that if data look roughly like following the assumed model, we will be better off than if they don't, so comparing the data with the assumed model makes some sense, but it has severe limitations. For example, if you take any well defined distance between distributions, it is not true in general that the LS-estimator will work the better, the closer the empirical residual distribution is to the normal distribution.

The task of model diagnostics is not to find out any "best" model for the errors, but rather to see whether the data have characteristics (such as outliers) that will lead to potentially misleading results, which also needs to take into account how in the real situation of interest results are going to be interpreted and used, and whether it for example is possible to monitor the data generating process in such a way that erroneous observations that are potential outliers can be reliably spotted and removed before analysis.