Linear regression with Cauchy distribution for errors

Question

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%

Answer 1

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)

Answer 2

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).