Regression when the OLS residuals are not normally distributed

Question

There are several threads on this site discussing how to determine if the OLS residuals are asymptotically normally distributed. Another way to evaluate the normality of the residuals with R code is provided in this excellent answer. This is another discussion on the practical difference between standardized and observed residuals.

But let's say the residuals are definitely not normally distributed, like in this example. Here we have several thousand observations and clearly we must reject the normally-distributed-residuals assumption. One way to address the problem is to employ some form of robust estimator as explained in the answer. However I am not limited to OLS and in facts I would like to understand the benefits of other glm or non-linear methodologies.

What is the most efficient way to model data violating the OLS normality of residuals assumption? Or at least what should be the first step to develop a sound regression analysis methodology?

Answer 1

The ordinary least squares estimate is still a reasonable estimator in the face of non-normal errors. In particular, the Gauss-Markov Theorem states that the ordinary least squares estimate is the best linear unbiased estimator (BLUE) of the regression coefficients ('Best' meaning optimal in terms of minimizing mean squared error)as long as the errors

(1) have mean zero

(2) are uncorrelated

(3) have constant variance

Notice there is no condition of normality here (or even any condition that the errors are IID).

The normality condition comes into play when you're trying to get confidence intervals and/or $p$-values. As @MichaelChernick mentions (+1, btw) you can use robust inference when the errors are non-normal as long as the departure from normality can be handled by the method - for example, (as we discussed in this thread) the Huber $M$-estimator can provide robust inference when the true error distribution is the mixture between normal and a long tailed distribution (which your example looks like) but may not be helpful for other departures from normality. One interesting possibility that Michael alludes to is bootstrapping to obtain confidence intervals for the OLS estimates and seeing how this compares with the Huber-based inference.

Edit: I often hear it said that you can rely on the Central Limit Theorem to take care of non-normal errors - this is not always true (I'm not just talking about counterexamples where the theorem fails). In the real data example the OP refers to, we have a large sample size but can see evidence of a long-tailed error distribution - in situations where you have long tailed errors, you can't necessarily rely on the Central Limit Theorem to give you approximately unbiased inference for realistic finite sample sizes. For example, if the errors follow a $t$-distribution with $2.01$ degrees of freedom (which is not clearly more long-tailed than the errors seen in the OP's data), the coefficient estimates are asymptotically normally distributed, but it takes much longer to "kick in" than it does for other shorter-tailed distributions.

Below, I demonstrate with a crude simulation in R that when $y_{i} = 1 + 2x_{i} + \varepsilon_i$, where $\varepsilon_{i} \sim t_{2.01}$, the sampling distribution of $\hat{\beta}_{1}$ is still quite long tailed even when the sample size is $n=4000$:

set.seed(5678)
B = matrix(0,1000,2)
for(i in 1:1000)
{
    x = rnorm(4000) 
    y = 1 + 2*x + rt(4000,2.01)
    g = lm(y~x)
    B[i,] = coef(g)
}
qqnorm(B[,2])
qqline(B[,2])

Answer 2

I think you want to look at all the properties of the residuals.

1. normality
2. constant variance
3. correlated to a covariate.
4. combinations of the above

If it is just 1 and it is due to heavytails or skewness due to one heavy tail, robust regression might be a good approach or possibly a transformation to normality. If it is a non-constant variance try a variance stabilizing transformation or attempt to model the variance function. If it is just 3 that suggests a different form of model involving that covariate. Whatever the problem bootstrapping the vectors or reiduals is always an option.

Answer 3

For non-normal conditions one would sometimes resort to robust regression, especially using the links to methods.

In order to present the context for non-normality it may help to review the assumptions for linear OLS regression, which are:

• Weak exogeneity. This essentially means that the predictor variables, x, can be treated as fixed values, rather than random variables. This means, for example, that the predictor variables are assumed to be error-free—that is, not contaminated with measurement errors. This assumption is the one that is most frequently violated and leads to errors as enumerated following this assumption list.

• Linearity. This means that the mean of the response variable is a linear combination of the parameters (regression coefficients) and the predictor variables. Note that this assumption is much less restrictive than it may at first seem. Because the predictor variables are treated as fixed values (see above), linearity is really only a restriction on the parameters. The predictor variables themselves can be arbitrarily transformed, and in fact multiple copies of the same underlying predictor variable can be added, each one transformed differently.

• Constant variance (a.k.a. homoscedasticity). This means that different values of the response variable have the same variance in their errors, regardless of the values of the predictor variables. In practice this assumption is invalid (i.e. the errors are heteroscedastic) if the response variable can vary over a wide scale. In order to check for heterogeneous error variance, or when a pattern of residuals violates model assumptions of homoscedasticity (error is equally variable around the 'best-fitting line' for all points of x), it is prudent to look for a "fanning effect" between residual error and predicted values. This is to say there will be a systematic change in the absolute or squared residuals when plotted against the predictive variables. Errors will not be evenly distributed across the regression line. Heteroscedasticity will result in the averaging over of distinguishable variances around the points to get a single variance that is inaccurately representing all the variances of the line. In effect, residuals appear clustered and spread apart on their predicted plots for larger and smaller values for points along the linear regression line, and the mean squared error for the model will be wrong.

• Independence of errors. This assumes that the errors of the response variables are uncorrelated with each other. (Actual statistical independence is a stronger condition than mere lack of correlation and is often not needed, although it can be exploited if it is known to hold. This latter can be examined with cluster analysis and correction for interaction.) Some methods (e.g. generalized least squares) are capable of handling correlated errors, although they typically require significantly more data unless some sort of regularization is used to bias the model towards assuming uncorrelated errors. Bayesian linear regression is a general way of handling this issue.

• The statistical relationship between the error terms and the regressors plays an important role in determining whether an estimation procedure has desirable sampling properties such as being unbiased and consistent.
  

• The arrangement, or probability distribution of the predictor variables x has a major influence on the precision of estimates of β. Sampling and design of experiments are highly developed subfields of statistics that provide guidance for collecting data in such a way to achieve a precise estimate of β.
  

As this answer illustrates, simulated Student's-$t$ distributed $y$-axis errors from a line lead to OLS regression lines with confidence intervals for slope and intercept that increase in size as the degrees of freedom ($df$) decrease. For $df=1$, Student's-$t$ is a Cauchy distribution and the confidence intervals for slope become $(-\infty,+\infty)$.

It is arbitrary to invoke the Cauchy distribution with respect to residuals in the sense that when the generating errors are Cauchy distributed, the OLS residuals from a spurious line through the data would be even less reliable, i.e., garbage in---garbage out. In those cases, one can use Theil-Sen regression regression. Theil-Sen is certainly more robust than OLS for non-normal residuals, e.g., Cauchy distributed error would not degrade the confidence intervals and unlike OLS is also a bivariate regression, however in the bivariate case it is still biased. Passing-Bablok regression can be more bivariate unbiased, but does not apply to negative regression slopes. It is most commonly used for methods comparison studies. One should mention Deming regression here, as unlike the Theil-Sen and Passing-Bablok regressions, it is an actual solution to the bivariate problem, but lacks the robustness of those other regressions. Robustness can be increased by truncating data to include the more central values, e.g., random sample consensus (RANSAC) is an iterative method to estimate parameters of a mathematical model from a set of observed data that contains outliers.

What then is bivariate regression? A lack of testing for bivariate nature of problems is the most frequent cause for OLS regression dilution and has been nicely presented elsewhere on this site. The concept of OLS bias in this context is not well recognized, see for example Frost and Thompson as presented by Longford et al. (2001), which refers the reader to other methods, expanding the regression model to acknowledge the variability in the $x$ variable, so that no bias arises$^1$. In other words, bivariate case regression sometimes cannot be ignored when both the $x$- and $y$-values are randomly distributed. The need for bivariate regression can be tested for by fitting an OLS regression line to the bijected residuals from an OLS regression of the data.

[Edit: Usually, residuals (the model minus the fitted function values) are plotted as versus the corresponding $x$-values. This will not show OLS bias. However, plotting residuals versus their sequence number: first, second, third, etc., that is, equidistantly, A.K.A as a bijected set, DOES show the bias]

Then, if the bijected OLS residuals have a non-zero slope, the problem is bivariate and the OLS regression of the data will have a slope magnitude that is too shallow, and an intercept that is too large in magnitude to be representative of the functional relationship between $x$ and $y$. In those cases, the least error linear estimator of $y$-values indeed would still be from OLS regression, and its R$^2$-value will be at a maximum possible value, but the OLS regression line will not represent the actual line function that relates the $x$ and $y$ random variables. As a counter example, when, as occurs among other problems in a time series with equidistant $x$-values, OLS of the raw data is not always inappropriate, it may represent the best $y=f(x)$ line, but is still subject to variable transformation, for example for count data, one would take the square root of the counts to convert the errors for Poisson distributed error to more normal conditions, and one should still check for non-zero slope of residuals.

1. Longford, N. T. (2001). "Correspondence". Journal of the Royal Statistical Society, Series A. 164: 565. doi:10.1111/1467-985x.00219

Answer 4

Macro (jsut above) stated the correct answer. Just some precision because I had the same question

The condition of normality of the residuals is useful when residuals are also homoskedastic. The result is then that OLS has the smallest variance between all of the estimator (linear OR non-linear).

The extended OLS assumptions:

1. $E(u|X_i = x) = 0$
2. $(X_i,Y_i), i=1,…,n,$ are i.i.d
3. Large outliers are rare
4. u is homoskedastic
5. u is distributed $N(0,σ^2)$

if 1-5 verified, then OLS has the smallest variance between all of the estimator (linear OR non-linear).

if only 1-4 verified, then by Gauss-Markov, OLS is the best linear (only !) estimator (BLUE).

Source : Stock and Watson , Econometrics + my course (EPFL, Econometrics)

Answer 5

My experience is completely in accord with Michael Chernick. Not only at times does applying a data transformation makes the modeling error normally distributed, it can also correct heteroskedasticity.

Sorry, but to suggest otherwise like gather an insane amount of data, or employ less efficient robust regression methods, is misguided, in my opinion, having practice this science/art.