# Linear regression with Cauchy distribution for errors

I have ran the below linear regression model and using the performance package in R I however checked whether the distribution of the residuals is normal. The performance package suggests I should be using a Cauchy distribution for the errors. In a search of stats.stackexchange and Google, it isn't clear how to do this. How can I model the data below using a Cauchy distribution for errors?

library(performance)

c(7L, 50L, 12L, 20L, 6L, 12L, 30L, 3L, 21L, 43L, 42L, 35L, 18L, 6L, 23L, 16L, 8L, 43L, 10L, 24L, 19L, 30L, 13L, 9L, 6L, 17L, 46L, 14L, 8L, 25L, 16L, 9L, 28L, 11L, 3L, 28L, 38L, 37L, 6L, 25L, 27L, 24L, 5L, 1L, 9L, 4L, 14L, 22L, 0L, 11L, 17L, 1L, 5L, 37L, 52L, 16L, 2L, 0L, 12L, 13L, 2L, 16L, 8L, 2L, 3L, 15L, 23L, 24L, 1L, 18L, 17L, 18L, 3L, 40L, 2L, 32L, 24L, 17L, 1L, 2L, 3L, 30L, 17L, 5L, 33L, 15L, 19L, 20L, 3L, 0L, 2L, 2L, 8L, 18L, 7L, 3L, 18L, 0L, 17L, 20L) -> dependent.var

c(4.66666666666667, 75, 28, 6, 1.83, 38.36, 80, 0, 14, 107, 137, 94.75, 36, 10.8666666666667, 44, 27, 32, 86, 52.8333333333333, 108, 76.5, 54, 26, 23.75, 11.75, 33.2133333333333, 100, 58, 50, 94, 32.25, 16, 33.75, 29.25, 7.75, 100, 98, 58.45, 4.58, 56, 59, 73.4166666666667, 6.16666666666667, 1, 53.79, 41.95, 43.25, 70.5, 0, 10, 3.25, 0, 14, 98, 112, 35, 0.25, 16.25, 30.83, 68, 1.25, 30.25, 13.25, 11.1, 1.5, 41, 45.17, 40, 6, 52.8566666666666, 43, 41, 3, 131, 0, 45.67, 74, 25.4166666666667, 0.25, 4.75, 14.58, 2.75694444444444, 32, 0, 92.25, 34, 66, 14, 1.75, 1.5, 1, 21.53, 4.08333333333333, 44.07, 55.9, 12, 20, 12.5, 48.1333333333333, 24.03) -> independent.var

lm(dependent.var ~ independent.var) -> model

check_distribution(model)

# # Distribution of Model Family
#
# Predicted Distribution of Residuals
#
#  Distribution Probability
#        cauchy         53%
#        normal         41%
#           chi          6%
#
# Predicted Distribution of Response
#
#                Distribution Probability
#  neg. binomial (zero-infl.)         59%
#               beta-binomial         34%
#                 half-cauchy          3%

• Does this package have any citations for the method used? My glance at it made it look like the answer is in the negative, and that leaves me reluctant to trust its output. The closest I found was the function documentation saying that it uses a random forest to classify the residuals into a distribution and that the algorithm for doing this might be improved in the future, but there was not a citation of a peer-reviewed article or conference paper, not even a working paper on arXiv. Anyone can upload a package to CRAN as long as it compiles; CRAN does not check the statistics aspects.
– Dave
Commented Apr 3, 2023 at 19:10
• What are these data? They seem to be only non-negative whole numbers. The normality of the residuals is generally the least important assumption in an OLS model. Even if the errors truly are Cauchy, the model fit should be fine. It's the SE's (& thus, CIs & p-values) that would be affected. There are robust methods to get a test, if that's what you want. What is the point of this model? Commented Apr 3, 2023 at 19:11
• Plotting the data makes me even more skeptical about the package output. The QQ plot is not great but is not so awful to suggest infinite variance, and a scatterplot of the data themselves looks like a fairly routine linear trend.
– Dave
Commented Apr 3, 2023 at 19:27
• @luciano that is not a good assumption to make for a general R package. as dave said, anyone can put code on CRAN, that doesn't mean the code implements a useful statistical method. Commented Apr 4, 2023 at 12:49
• @wzbillings I completely agree by the way with your assessment that it's not good idea to use the check_distribution function blindly (at all, probably). So not sure what point I'm trying to make other than its authors are not random people who don't know what they are doing. It was reminded of fortune #386: "If we put in a function into rstan that dropped chains, people would use it.". Commented Apr 4, 2023 at 17:44

This question raises a couple of interesting points:

• In a comment, @Dave points out that the normal QQ plot is "not so awful to suggest infinite variance".
• In an answer (+1), @wzbillings points out that it's not clear what theory (if any) the performance package uses to justify the Cauchy distribution for the errors.

As @Dave advises, let's look at the data first. (I've renamed the independent variable $$x$$ and the dependent variable $$y$$ for clarity.)

Clearly both variables are nonnegative (and in fact the values are mostly integers but the OP hasn't provided many details about the data). The relationship between $$x$$ and $$y$$ is reasonably linear and the variability of $$y$$ doesn't appear to increase for large $$x$$. Except near the origin both variables are constrained and that "squeezes" the residuals: their variance is lower for $$x$$ close to 0 and as a result the distribution of the residuals is heavy tailed.

Rather than choose between the Normal and the Cauchy distribution for the errors, let's fit a simple linear regression with $$t$$-distributed errors. And we will use theory to explore what degrees of freedom, $$\nu$$, are most consistent with the data. (Recall that the Cauchy is equivalent to $$t$$ distribution with $$\nu = 1$$ degree of freedom and the Normal is equivalent to $$t$$ with infinitely many degrees of freedom.)

The idea is to compute (numerically) the profile likelihood $$L(\nu)$$ for the degrees of freedom parameter $$\nu$$.

We start by writing down the likelihood of all four parameters in a simple linear regression with $$t$$-distributed errors:

\begin{aligned} L(\beta_0,\beta_1,\sigma,\nu) = \prod_i f_\nu\left(\frac{y - (\beta_0 + \beta_1x)}{\sigma}\right) \times \frac{1}{\sigma} \end{aligned} where $$f_\nu$$ is the density of the $$t$$ distribution with $$\nu$$ degrees of freedom and $$1/\sigma$$ is the Jacobian of the transformation $$y_0 = (y - \mu)/\sigma = (y - \beta_0 - \beta_1x)/\sigma$$.

Then we find the profile likelihood $$L(\nu)$$ of the degrees of freedom parameter by maximimizing the likelihood $$L(\beta_0,\beta_1,\sigma,\nu)$$ with respect to $$\beta_0,\beta_1$$ and $$\sigma$$. (Actually, we minimize the negative log likelihood; see R code below.)

\begin{aligned} L(\nu) = \max_{\beta_0,\beta_1,\sigma} L(\beta_0,\beta_1,\sigma,\nu) \end{aligned}

And here is the plot of the profile likelihood.

A couple of conclusions from this analysis:

• The profile likelihood $$L(\nu)$$ is maximized at $$\nu = 4$$ but anything between 3 and 5 degrees of freedom is a good fit to the data.
• As @Dave claimed, the Normal distribution ($$\nu \approx 30$$) is a better fit than the Cauchy distribution ($$\nu = 1$$).

References

[1] In All Likelihood: Statistical Modelling And Inference Using Likelihood. Y. Pawitan. Oxford University Press (2013)
[2] A note by @kjetilbhalvorsen about the dangers of maximum likelihood estimation for the degrees of freedom of the $$t$$ distribution. The warning is about small sample sizes in particular while in this case $$n = 100$$, so the profile likelihood approach is probably okay.

# Requires data x and y
n <- length(y)

m0 <- lm(y ~ x)

library("nloptr")

minimize <- function(x0, func, lb = NULL, ub = NULL) {
opts <- list(
"algorithm" = "NLOPT_LN_SBPLX",
"xtol_rel" = 1.0e-6
)
a <- nloptr(x0, func, lb = lb, ub = ub, opts = opts)
list(x = a$$solution, fx = a$$objective)
}

negloglik <- function(beta0, beta1, sigma, nu) {
mu <- beta0 + beta1 * x
y0 <- (y - mu) / sigma
# Calculate the log likelihood for y0
ll <- sum(dt(y0, nu, log = TRUE))
# Add the log of the Jacobian
ll <- ll - n * log(sigma)
# Calculate the negative log likelihood foy y
-ll
}

beta0.hat <- m0$$coef[1] beta1.hat <- m0$$coef[2]
sigma.hat <- sigma(m0)

profile <- function(nu) {
sapply(nu, function(nu) {
soln <- minimize(
c(beta0.hat, beta1.hat, sigma.hat),
function(params) negloglik(params[1], params[2], params[3], nu),
lb = c(-Inf, -Inf, 0),
ub = c(Inf, Inf, Inf)
)
soln\$fx
})
}

nu <- seq(1, 30, by = 0.25)
nll <- profile(nu)
ll <- exp(min(nll) - nll)

plot(nu, ll,
type = "l",
xlab = quote(paste("degrees of freedom, ", nu)),
ylab = quote(L(nu)),
main =  quote(paste("Profile likelihood L(", nu, ")"))
)
abline(h = 0.15, lwd = 0.3)


I don't really like this approach to determining the model family, I think it should be set a priori instead of using this machine learning method that this package uses with no citations to methodological papers. But that isn't the question.

You could use the function heavy::heavyLm() with a t-distribution family and manually set the degrees of freedom to 1. (Note that this package is not on CRAN: https://github.com/faosorios/heavy)

Update: apparently the package also has a built-in Cauchy family for heavyLm(), so you could use that also. They should work out the same (to a meaningful level of precision anyways).

• Stack Overflow supports your suggestion to use the heavy package.
– Dave
Commented Apr 3, 2023 at 19:27
• The criticism about the performance::check_distribution function is not mis-misplaces but still the documentation states: "This function is somewhat experimental and might be improved in future releases. The final decision on the model-family should also be based on theoretical aspects and other information about the data and the model." So fair to say that the package authors are clear about the pitfalls. Commented Apr 4, 2023 at 9:38
• @dipetkov I'm not sure I agree, I think that "this is experimental and may be improved" and "this method has no empirical or theoretical evidence to suggest it does what it is supposed" are not the same thing Commented Apr 4, 2023 at 12:48