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 residuals from an OLS regression of the data. Then, if the 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.
- Longford, N. T. (2001). "Correspondence". Journal of the Royal Statistical Society, Series A. 164: 565. doi:10.1111/1467-985x.00219